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1.  OVERVIEW 

This  document  is  the  finsl  report  for  Fusion  of  Correlated  Sensor  Observations,  an  effort 
cffluducte^*  by  igaman  Sdences  Corporaticm  (KSC)  under  Rome  Laboratory’s  Broad  Agency 
Announcement  (BAA)  number  90-04.  This  effort  was  performed  for 

Dr.  James  H.  Michels 
Rome  Ldxsratory 
RWOCTM 

Griffiss  AFB,  NY  13441 
(315)  330-4432 


1.1  BtMdtground 

Over  the  last  few  years,  Dr.  James  Michels  of  RL/OCTM  has  been  investigating  an 
(viginal  approach  in  signal  processing.  He  has  described  this  approach  and  results  to  date  in  a 
series  of  Rome  Laboratory  Technical  Reports  (Michels  89,  90a,  90b,  91,  92a,  and  92b).  This 
investigation  has  developed: 

•  A  synthesis  procedure  for  simulating  multichannel  autoregressive  (AR)  processes  in 
which  intertemporal  and  interchannel  ctxrelations  are  controlled  parametrically. 

•  Extensioos  to  this  synthesis  procedure  to  handle  non-Gaussian  Spherically  Invariant 
Random  Process  (SJRJP)  noise  and  unconstrained  quadrature  cmnponents. 

•  Statistical  diagnostics  for  analyzing  the  expected  performance  of  various  multichaimel 
estimating  algorithms. 

•  A  multichannel  signal  detection  algorithm  based  cm  a  generalized  likelihood  ratio  using 
an  innovations  iq^oach.  Two  rstios  currently  exist,  one  for  Gaussian  processes  and  one 
for  non-Gaussian  SIRPs. 

1.1.1  The  Muttkhannel  Signal  Processing  Simulation  System 

Kaman  Sciences  has  designed  and  implemented  a  Graphical  User  Interface-based  fourth 
generation  system  fw  analyzing  the  algorithms  developed  during  this  investigation.  This 
Multichannel  Signal  Processing  Simulation  System  (MSPSS)  is  described  in  the  Kaman 
Sciences  report  Man  Machine  Interface  Experiment  (Kaman  91a). 

The  MSPSS  is  com|xised  of  two  major  subsystems,  a  menu-based  subsystem  and  the 
User  Front-end  Interface  (UFI)  based  subsy^m.  The  menu-based  subsystmn  interacts  with  the 
user  by  a  series  of  menus  and  pronqns  diat  can  be  displayed  on  a  line-oriented  "glass 
terminal."  It  is  implemented  as  a  collecdon  of  Fortran  programs  providing  the  desired  signal 
{vocessing  capabilities.  Multichannel  data  is  passed  between  diese  programs  by  user-defined 
files.  The  menus  themselves  are  unix  c-shells  where  the  bottom-level  menu  automatically 
compiles.  Units,  and  executes  the  desired  program.  Figure  1-1  shows  an  example  of  the 
invocatitm  of  a  Fortran  program  from  the  menu-based  subsystem. 

The  structure  of  the  menu-based  subsystem  provides  the  user  with  a  great  deal  of 
analyticai  flexibiUty.  The  individual  Fortran  programs  provide  certain  analytical  capabiUties, 
but  minimal  constraints  are  imposed  on  the  order  in  which  the  user  can  invoke  these 
functions.  A  number  of  diagnostic  capabiUties  are  provided,  including  parameter  estimation 
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The  user  invotes  the  menu-based 

%  multichannel  system  from  the  Unix  shell 

Mokj-Oumnel  Detecdoa  Algorithm 

MAIN  MENU 

Sii^  Mhhirinniiel 

Ml  -  Rocess  Synibesis  MU  -  Ptaoett  Syatfaem 
M2  -  Hlteriag  Metbob  M12  -  Hlteriqg  Methods 
M3  -  Diagnostics  M13  -  Diagnostics 

LI  -  List  data  files  L2  -  show  the  change  log 
Enter  a  «»"«"««<<  or  Q  to  quit:  MU  The  user  chooses  a  submenu 

MULTI-CHANNEL  (M/O  PROCESS  SYNTHESIS  MENU 

MO  -  Generate  Gaussian  noise 

Ml  -  Method  1  (MCI)  process  synthesis 

M2  -  Method  2  (MC2)  process  synthesis 

M3  -  Method  2  process  synthesis  widi  unconstrained  quadrature  components 

M4  -  State  qtaoe  process  synthesis 

MS  -  Apply  Levinson-Wiggins-RotMnsoo  algorithm 

M6  -  Display  Multi-ChaiBiel  (M/Q  signal  stats 

M7  -PlotM/Cdaia 

M8  -  Feifonn  Nnttall-Strand  or  ^^dn-Morf  estimmion 

M9  -  Perform  Yule-WalkBr  estimation 

MIO  -  EstimalB  AR  psameteis  for  each  channel 

MU  -  EstimalB  covariance 

M12  -  Perform  estimatioa  of  stale  space  patammers 

M13  —  Perform  l/UC  correlation  (lerqxiral) 

M14  -  Perform  M/C  cooelation  (ensemUe) 

MIS  -  Perform  MiC  correlatioa  (quadrature,  tBnq)oral) 

M16  -  Perform  M/C  correlation  (quadrature,  ensemble) 

M17  -  Plot  ergodic  series  equation 

M18  -  Split  a  M/C  irqnt  M24  -  Join  iigaits 

M19  -  Add  M/C  signals  M2S  -  Subtract  M/C  signals 

M20  *>  Convert  data  lo  ASCH  file  M26*~  Convert  ASCII  file  to  data 

M21  -  Convert  AR  to  stale  space 

M22  -  Diqriay  oomplea  data  M27  -  Display  coefficient  file 

M23  -  Catenate  M/C  signals  M28  -  Sum  channels 

M29  -  Perform  SVD 

LI  —  List  data  files  L2  -  dbaw  the  change  log 
*  denotes  options  that  are  not  yet  avsHsUe 
Enter  a  command  or  Q  to  lenn  to  the  MAIN  menu:  M4 
Tjniring  program  mcssl,  please  stand  by  A  Fortran  program  is  called 

f77  -O  -C  -pipe  -cg89  sund/incssl.o  -Lsut^  -Itnc  -Ivec  -Imat  -Imc 
•lihaniefoiarlii)/lam/lib/san4  -Ifol  -s  -o  bin/sun4foicssl 
Synthesires  stale  qrace  process 
Version  1.4 

Do  yon  want  to  set  the  pseudo  random  generator  seeds  (y/n)  [N]  ?  _ 

The  user  interacts  with  the  Fortran  program 


Figure  1*1:  Invoking  a  Fortran  Program  from  the  Menu-Based  System 
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for  certain  models,  correlation  function  estimation,  and  graphical  display  of  data.  Thus  the 
user  is  {»'Ovided  with  a  system  that  supports  an  exploratory  style  of  analysis. 

The  UFI-based  subsystem  is  designed  for  more  repetitive  analyses  characteristic  of  the 
determination  of  detection  probabilities  for  given  false  alarm  probabilities  and  given 
algorithms.  In  such  analyses,  the  user  needs  to  invoke  certain  signal  processing  synthesis  and 
analysis  funcdtms  in  a  {xredetermined  ordff  with  certain  parameters.  The  UFI-based 
subsystem  si^Tpcrts  this  need  by  allowing  the  user  to  construct  an  ''e;q)eriment"  specified  by 
various  algoridims  and  parameters.  Once  an  expaiment  has  been  completely  described,  a 
computer  program  for  p^orming  an  experiment  is  automatically  generated,  compiled,  and 
executed.  Since  these  are  simulation  experiments,  the  generated  program  can  run  for  quite 
some  time.  Typically,  a  single  experiment  will  be  used  to  analyze  the  performance  of  the 
signal  detection  algorithm  in  terms  of  the  false  alarm  probability  and  the  probability  of 
detecticm.  A  detailed  example  of  the  use  of  this  subsystem  is  provided  in  the  Software  User’s 
Manual  for  the  Multichannel  Signal  Processing  Simulation  System  (Kaman  92a). 

1.12  Recent  Research 

Some  research  performed  in  parallel  with  the  development  and  use  of  the  MSPSS  has 
some  interesting  implications  for  Dr.  Michels’  approach. 

RL  has  been  sponsoring  work  in  sensor  fusion  (Al-Ibrahim  91,  Chair  86,  and  Hoballah 
89).  In  the  sensor  fusion  problem  several  senstxs,  perhaps  geographically  distributed,  each 
process  a  signal  they  receive.  These  sensors  send  their  outputs  to  a  centrally  located  fusion 
center  which  dien  must  decide  if  a  sigiutl  is  present.  The  fusion  center  combines  the  results  of 
each  of  the  individual  sensors  to  form  a  decision  for  the  overall  system.  One  can  view  each 
channel  in  a  multichannel  processing  problem  as  catrespaadiag  to  the  input  to  a  sensor  in  the 
sensor  fusion  problem.  Thus.  Dr.  Michels’  algoridun  is  a  promising  aj^oach  for  the  sensor 
fusion  problmn,  but  further  research  is  needed  to  determine  quantitative  benefits. 

RL  has  been  sponsoring  some  work  on  Kalman  filtering  (Roman  93).  In  Dr.  Michels’ 
software  detection  algorithm,  linear  filters  are  ai^lied  to  the  multichannel  radar  returns  before 
calculating  the  loglikelihood  ratio  upon  which  the  signal  detection  is  based.  The  tapped  delay 
line  filters  in  use  when  this  effort  began  can  be  replaced  by  other  filter  structures.  In 
particular,  a  Kalman  filter  is  iq>plicable.  Once  again  additional  research  is  needed  to  determine 
the  benefits  of  using  the  new  Kalman  filter  bdng  sptmsored  by  RL  in  Dr.  Michels’  signal 
detecticm  algorithm. 

Signal  detection  algcxithms  are  traditionally  evaluated  in  terms  of  the  false  alarm 
probability  (the  percentage  of  instances  in  which  the  algorithm  results  in  a  decision  that  a 
signal  is  peseat  when  no  signal  is  in  fact  present)  and  the  probability  of  detection  (the 
probability  that  the  algcxithm  will  register  a  signal  when  in  fact  a  sigmd  is  presort).  In  a 
simulaticm  app’oach  for  evaluating  a  signal  detection  algorithm,  a  very  large  number  of 
samples  needs  to  be  used  to  accurately  explore  the  tails  of  probability  distributions,  since  the 
false  alarm  probability  is  itself  a  "tsul  jx^obability.”  Recent  research  has  proposed 
approximating  these  tails  based  on  a  Pareto  distribution  (Chakravarthi  92  and  Rangaswamy 
93).  This  qrproach  promises  at  least  an  order  of  magnitude  improvement  in  throughput  in  Dr. 
Michels’  studies,  but,  once  again,  further  research  is  needed  to  validate  this  claim. 

Signal,  dutmr,  and  noise  are  often  modeled  based  on  certain  probability  distributions.  A 
powerful  algorithm  has  recently  been  developed  for  determining  the  probability  distribution 
from  which  a  small  random  sample  is  drawn  (Ozturk  90a,  90b.  90c).  The  Ozturk  algorithm 
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provides  the  MSPSS  with  an  imporant  additional  diagnostic  capability. 

To  explore  the  implications  of  sensor  fusion,  Kalman  filtering,  the  Pareto  distribution 
af^jroech,  and  the  Ozturk  algorithm  for  Dr.  Michels’  signal  detectitm  algorithm,  new 
capabilities  were  added  to  the  MSPSS.  The  result  of  incorporating  these  cl^}abilities  was  an 
«>nhaiiraH  simulstitm  System  for  exploring  the  performance  of  a  sensor  fusion  approach  using 
an  ionovatLtms-based  t^dinique. 

IJt  Objective 

The  ob^ctive  of  tins  effort  was  to  incorporate  additional  functionality  into  the 
Multichannel  Signal  Processing  Simulation  System  (MSPSS)  that  Kaman  Sciences 
Corparation  has  developed  to  support  the  analysis  described  in  Multichannel  Detection  Using 
the  Discrete-Time  Mo^l-Based  Innovations  Approach,  (Michels  91).  Additional  fimetions 
included  the  following: 

•  The  ability  to  support  analyse,  of  a  sensor  fusion  algorithm 

•  A  multichannel  Kalman  filter 

•  The  use  of  extreme  value  thetx^  to  set  thresholds  accurately  with  much  lower  sample 

sizes 

•  The  Ozturk  algtxithm  for  determining  a  probability  distribution  from  a  random  sample. 

U  Tasks  Performed 

This  effort  consisted  of  four  major  tasks.  Some  related  minor  additional  tasks  were  also 
performed  under  this  effort 

1J.1  Task  1:  Sensor  Fusion  Capabilities 

Under  tius  task,  we  added  two  programs  to  the  menu-based  subsystem  and  four  analysis 
sequences  to  the  UFI-based  subsystem.  One  of  the  menu-based  programs  estimates  the 
covaiiance  matrix,  and  the  other  estimates  Autoregressive  coefficients  for  each  channel.  Two 
of  tile  UFI-based  analysis  sequences  test  the  Autoregressive  (AR)  confidents  and  covariance 
estimatiem  algorithms.  The  other  two  sequences  analyze  sensor  fusion  detection  algorithms; 
they  differ  depending  upon  whether  temporally  correlated  noise  (or  clutter)  is  assumed  to  be 
present  or  not 

These  new  modules,  in  combination  with  previously  existing  modules,  provide  the  user 
with  the  capability  to  analyze  single  channel  processes  within  each  element  of  a  vector 
process.  This  capability  allows  the  user  to 

•  Form  scalar  linear  estimates  of  AR  coeffidents  for  each  channel  process 

•  Subsequently  process  each  single  channel  innovations  process  at  a  fusion  center  to 

remove  residual  correlations  across  channels. 

U.2  Task  2:  Kalman  Filtering 

When  this  effort  began,  a  software  filter  and  algorithm  implementation  of  a  multichannel 
Kalman  filter  was  being  developed  for  Rome  Labwatory  under  contract  F30602'92-C-0081, 
titled  State-Space  Models  for  Multichannel  Detection  (Roman  93).  This  task  implemented  in 
tim  MSPSS  the  algtxithms  developed  under  that  contract.  Kaman  Sciences  began  by 
examining  the  source  code  and  the  reports  produced  by  the  RL  contract  developing  the 
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multichannel  Kalman  filter.  We  then  developed  code  confmming  to  our  conventions  and 
interfaces  based  on  the  algorithms  developed  in  the  RL  contract. 

These  algorithms  were  first  implemented  in  the  menu-based  subsystem  of  the  MSPSS. 
They  were  then  implemented  in  the  UFI-based  subsystem.  New  capabilities  consisted  of 

•  The  synthesis  of  state  space  processes  in  the  menu-based  system 

•  Tim  estimation  of  state  space  model  parameters 

•  The  implementation  of  a  Kalman  filter 

•  The  implementation  of  a  signal  detection  algorithm  based  on  the  state  space  model  and 

Kalman  filter. 

1  Task  3:  Extreme  Value  Method 

Under  this  task,  the  extreme  value  theory  method  for  the  threshold  level  selection, 
described  in  (Rangaswamy  93),  was  implemented.  This  method  is  included  as  an  alternative 
to  the  previous  threshold  selection  method  contained  in  the  detection  modules  of  MSPSS 
software. 

Kaman  Sciences  first  developed  software  implementing  the  extreme  value  theory  as  a 
module  in  the  diagnostic  menu-based  software.  This  only  involved  adding  an  additional 
likelihood  function  module,  since  the  previous  likelihood  function  module  calculates  threshold 
values  and  detection  probabilities  as  well.  Once  it  was  tested,  the  extreme  value  module  was 
then  modified  to  run  under  the  User  Front-&id  Interface  (UFI).  Kaman  Sciences  wrote 
additional  detection  sequences  that  parallel  the  previous  functionality,  but  use  the  extreme 
value  method. 

U.4  Task  4:  Ozturk  Algorithm 

Ozturk’s  algwithm,  as  described  in  (Ozturk  90a,  90b,  90c),  (Shah  93),  and  (Slaski  93), 
was  implemented  under  this  effort.  This  implementation  resulted  in  two  programs  under  the 
diagnostic  menu-based  system.  Ozturk’s  algorithm  is  ^plicable  to  single-channel  real  data. 
The  smaller  program  provides  a  firont-end  fcx*  converting  many  realizations  of  multichannel 
data  into  a  single  realization  of  single  channel  real  data  by  calculating  a  quadratic  form. 

The  second  program  calculates  Ozturk’s  statistic.  It  performs  both  a  goodness-of-fit  test 
and  estimates  a  probability  distribution  from  the  input  data.  The  distribution  of  Ozturk’s 
statistic  is  not  known  in  closed  form.  Thus,  this  program  includes  an  extensive  Monte-Carlo 
c^>ability  for  determining  the  expected  value  and  confidence  regions  fm*  the  statistic  under  a 
wide  range  of  null  hypotheses.  An  important  capability  of  Ozturk’s  algtxithm  is  the  ability  to 
graphically  display  the  results.  The  ^goithm  combines  a  formal  statistical  test  with  more 
intuitive  griq}hical  analysis.  Our  implementation  includes  an  integrated  grt^hical  display. 
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2.  THE  MULTICHANNEL  SIGNAL  PROCESSING  SIMULATION 
SYSTEM 

Tbe  Multichannel  Signal  Processing  Simulation  System  (MSPSS)  was  developed  by 
Kaman  Sciences  to  sui^)ort  certain  RL/OCTM  research.  This  section  briefly  overviews  the 
problems  that  can  be  analyzed  with  the  MSPSS,  the  software  design,  and  how  new 
capabilities  are  added  to  the  system. 

2.1  Signal  Detection  as  Statistical  Hypothesis  Testing 

The  MSPSS  simulates  correlated  random  vectors  which  characterize  several  channels  of 
data.  Hw  data  gathered  at  any  point  in  time  consists  of  a  real,  imagiiuuy,  or  complex  vector 
witii  a  component  fn*  each  channel  A  single  realization  or  trial  of  the  stochastic  process 
generating  tto  data  consists  of  an  wdered  sequence  of  vectms. 

The  signal  detection  problem  can  be  cast  in  the  form  of  a  problem  in  statistical 
hypothesis  testing.  Let  x  denote  a  vecto-  stochastic  process  representing  the  radar  return,  s 
denote  a  signal,  c  denote  the  clutter,  and  n  denote  white  noise.  Consider  deciding  between  the 
null  hypotiiesis 

•  Ho:  X  me  +n 
and  the  alternative 

•  Hi:  X  ms  +  C  n. 

To  accq)t  the  alternative  hypothesis  is  to  decide  that  a  signal  is  present. 

Neyman-Pearson  theory  is  generally  used  to  evaluate  a  decision  algorithm  in  statistical 
hypothesis  testing.  A  test  statistic  is  defined  for  summarizing  the  data  in  a  single  realization  of 
tte  stochastic  process.  If  this  statistic  is  in  some  critical  tegimi,  then  the  alternative  hypothesis 
is  aocq>ted.  Otiierwise  the  null  hypothesis  is  accq>ted.  The  probability  of  accq)ting  the 
alternative  when  tiie  null  is  true  is  called  tl»  significance  level;  this  i»robability  should  be  as 
low  as  possible.  The  probability  of  deciding  in  flavor  of  the  alternative  when  in  fact  it  is  true 
is  known  as  the  power  of  the  test.  Statistics  are  designed  to  maximize  the  power  of  the  test 
for  a  given  significance  level  and  sample  size.  Alternatively,  they  are  designed  to  minimize 
the  sanq}le  size  for  a  given  significance  level  and  power. 

Because  of  the  physical  interpretation  of  this  particular  statistical  test,  different 
terminology  is  used  in  radar  tiieory.  The  probability  of  erroneously  concluding  a  signal  is 
present  is  the  false  alarm  probab^ty.  So  the  false  alarm  {vobability  is  the  same  as  the 
significance  level.  The  probability  of  detecting  a  signal  wlmn  tiiete  is  in  fact  a  signal  is  the 
probability  of  detection,  or  the  detection  probability.  Therefore  the  detection  probability  and 
tile  power  of  the  test  are  alternate  names  for  the  same  jx’obability. 

12  Hie  Signal  Detection  Algorithm 

Figure  2-1  shows  the  general  structure  of  the  signal  detection  algorithm  implemented  in 
the  MSPSS.  This  is  a  model-based  iimovations  approach  to  signal  detection,  llie  filters  are 
derived  from  certain  models  of  the  signal,  clutter,  and  noise.  If  the  model  for  the  null 
hypotiiesis  filter  accurately  describes  the  sum  of  clutter  and  noise,  and  the  input  is  clutter  and 
noise,  tiien  the  ouqiut  of  the  null  hypothesis  filter,  yo(n ),  will  be  uncorrelated  both  in  time 
and  across  channels.  The  outputs  are  estimates  of  the  innovations  for  the  model  underlying  the 
filter.  Similarly,  if  the  model  for  the  alternative  hypothesis  filter  accurately  describes  the  sum 
of  the  signal,  clutter,  and  noise,  the  output  y,  (n )  will  be  an  estimate  of  the  innovations  for  the 
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altem«tive  hypothesis  model,  uncorrelated  both  in  time  and  aaoss  channels,  when  the  input  to 
the  filter  contains  signal,  clutter,  and  noise. 


Figure  2-1:  System  Structure  of  the  Signal  Detection  Algorithm 


The  loglikelihood  statistic  is  based  oa  the  estimates  of  the  innovadoas  produt^  by  the 
two  filters,  ff  this  statistic  is  above  a  jnedefined  threshold,  one  deades  a  signal  is  ;»esent. 
Otherwise,  the  alternative  hypothesis  is  rejected,  and  (me  decides  no  signal  is  jnesent 

The  structure  of  this  algorithm  is  very  general  and  raises  many  questions.  Different 
can  be  used  for  the  signal,  clutter,  and  noise.  Some  of  these  models  can  be  used  to 
approximate  other  models.  For  example,  if  the  signal  and  clutter  are  Autoregressiw  (W 
processes  and  the  noise  is  white,  their  sum  is  an  Autoregressive  Moving  Average  (ARMA; 
process.  An  ARMA  process  can  be  approximated  by  a  high  order  AR  process.  the  other 
a  State  Space  model  can  exactly  describe  an  ARMA  process.  How  does  this  choice  of 
model  impact  the  performance  of  the  algorithm? 

Models  may  also  be  chosen  on  the  basis  of  processing  considerations.  For  example, 
suppose  the  sum  of  the  signal,  clutter,  and  noise  are  approximated  by  a  vector  AR  process 
dneribed  by  Equation  2-1: 

x{n)»-  a"  {k)x{n -k)  +  t(n),  (2-1) 

k  m  I 

where  x  ( 1 ),  x  (2) . x  (^ )  is  the  AR  process  and  e(n )  is  a  random  (possibly  complex) 

mean  vector  with  covariance  matrix  £.  If  the  off-diagonal  elements  of  £  arc  non-zero,  the  time 
x(rt)  are  corelated  across  channels.  The  matrix  AR  coefficients  A"  (k)  lead  to 
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correlatiw  of  the  i»-ocess  both  in  time  and  across  channels.  The  off-diagonal  components  of 
the  matrix  AR  coeffiicents  may  be  non-zero.  If  they  were  all  zero.  Equation  2-1  would  be 
equivalent  to  Equadoa  2-2; 


Xjin) 


"  L  aj(t)x(n  -k)-t-ej(n). 

tml 


(2-2) 


j  s  1,  2,  ....  /  (/  is  die  number  of  channels).  The  model  given  by  Equaticm  2-2  has  the 
advantage  dut  the  catreqxmding  tapped  delay  line  filters  Fo  and  Fi  can  process  each  channel 
separately  in  die  signal  detection  algorithm.  This  model  can  be  implemented  with  considerable 
parallelism,  thus  gaining  an  increase  in  processing  speed.  But  is  this  increase  in  processing 
speed  accompanied  by  a  decrease  in  the  detection  probability  for  given  false  alarm 
probabilities? 

In  fxactice,  model  parameters  will  not  be  known  exacdy.  This  raises  a  host  of  new 
questims.  For  each  model,  what  parameter  estimation  algorithms  are  available?  How  do  they 
perfonn  in  various  circumstances?  For  example.  Dr.  Michels  has  found  that  the  conditions  that 
increase  the  accuracy  of  estimates  of  the  AR  wei^ts  increase  the  variability  of  estimates  of 
the  driving  noise  covariance  matrix  (Michels  92b).  How  robust  is  the  signal  detection 
algorithm  to  errors  in  parameter  estimation  and  model  mismatch?  These  questions  probably  do 
not  have  one  simple  answer.  Rather,  the  signal  detection  performance  of  any  one  instantiation 
of  the  algorithm  shown  in  Figure  2-1  will  vary  with  process  characteristics  such  as  temporal 
correlation.  A  full  analysis  of  the  model-based  innovatitxis  approach  to  signal  detection  will 
have  to  e3q>lore  the  variability  in  performance  with  changes  in  the  underlying  stochastic 
processes.  Detection  probabilities  and  duesholds  will  need  to  be  known  for  given  false  alarm 
probabilities  and  sample  sizes  in  implementing  the  signal  detection  algorithm.  A  complete 
analysis  of  the  signal  detection  algoritisym  will  also  provide  these  needed  parameters. 

The  Multichannel  Signal  Processing  Simulatitm  System  provides  the  analyst  with  tools  to 
answer  tiiese  questions  via  Monte-Carlo  simulatitm.  Capabilities  are  provided  to  synthesize 
signal,  clutter,  noise,  and  their  stun  as  des^ibed  by  a  variety  of  models.  The  user  specifies  the 
number  of  jMrocess  realizations  to  generate,  as  well  as  the  number  of  time  samples  in  each 
realization.  Diagnostic  capabilities  are  provided  for  analyzing  the  realizatitms  of  the  generated 
stochastic  imocesses.  The  user  can  calculate  means,  variances,  covariance  matrix  estimates, 
•  correlation  functions,  model  parameter  estimates,  and  Fast  Fourier  Transformations  (FFTs). 
Grq)hical  aq)abilities  are  provided  for  users  operating  under  a  X  Windows  interface. 
Functions  are  provided  for  manipulating  radar  returns  such  as  adding  them,  splitting  a 
multichannel  return  into  several  single  channel  returns,  combining  single  channel  returns  into 
multichannel  returns,  and  so  tm. 

Ifinally,  die  MSPSS  implements  the  signal  detection  algorithm  described  earlier.  Trqiped 
delay  line  and  Kalman  filters  are  provided,  with  user-defiiutble  parameters,  for  removing 
temporal  and  cross-channel  correlation.  Several  loglikelihood  statistics  can  be  calculated  for 
different  models.  Thresholds,  detection  probabilities,  and  variances  in  estimates  of  these 
parameters  can  be  calculated. 


23  Using  the  MSPSS 

As  was  explained  in  Section  1.1.1,  two  major  subsytems  comprise  the  Multchannel 
Signal  Ftocessing  Simulation  System  (MSPSS),  the  menu-based  subsystem  and  the  UFI-based 
subsystem.  The  use  of  the  UH-based  subsystem  is  extensively  described  in  the  Software 
User’s  Manual  for  the  Multichannel  Signal  Processing  Simulation  System  (Kaman  92a)  and 
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will  not  be  further  examined  here. 

As  shown  in  Figure  1-1,  invoking  the  menu-based  subsystem  presents  the  user  with  a 
choice  of  two  sets  of  three  menus.  These  menus  provide  the  user  with  a  wide  range  of 
capabilities.  The  bottom-level  menu  initiates  execudmi  of  a  Fortran  program.  The  programs 
in  this  system  share  common  conventions  for  user  into'action.  Figures  2-2  and  2-3  show  an 
example.  The  uso*  is  prompted  for  various  parameters  specifying  the  synthesis  or  analysis 
{vocess.  In  the  example  the  user  describes  tte  synthesis  of  an  Autor^ressive  ix’ocess  by  a 
"shaping  function"  ^jproach.  The  user  can  receive  help  messages  for  prompts  that  are  not 
self-evident,  but  a  ^  knowledge  of  the  processing  depends  on  an  understanding  of  the 
underlying  mathematical  models. 

The  menu-based  subsystem  provides  the  user  with  a  graphical  display  capability,  as 
illustrated  in  Figures  2-4  throu^  2-6.  This  capability  requires  a  X  Windows  interface.  The 
plotting  rmitine  first  asks  die  user  for  the  source  of  Ae  input  file  and  various  parameters  fcx 
the  graph  (Figure  2-4).  An  user-friendly  windowing  interface  then  appears,  permitting  the  user 
to  manipulate  the  plot  (Figure  2-5).  Outputs  sent  to  the  screen  can  also  be  saved  as  a 
Postscript  file  or  printed  (Figure  2-6). 

An  oudine  of  the  steps  required  to  add  capabilities  to  the  MSPSS  is  provided  in 
(Vieimeau  93).  Both  the  menu-based  and  UFI-based  subsystems  are  described. 
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The  menus  are  displayed 


Enter  a  oommaDd  or  Q  to  return  to  the  MAIN  meim:  ai2 
Linldiig  progrm  msirp.  please  stan.^  by 
‘biii/niii4Aiiatip'  is  iq>  to  date. 

MSIRP  -  MabkhaDiiel  AR  Process  Geoeratiao  with  SIRPs 
Venion  1.7 

Do  you  want  to  set  tbe  pseudo  random  generator  seeds  (y/n)  [N]  ? 

Number  d  rJiMinaU  ?  A  Carriage  Return  obtains  help 

Enter  the  number  of  channels  to  be  genemied  in  eadi  output  signal  or 
enter  a  zero  to  quit  die  program. 

Number  of  channels  ?  1 

Number  of  points  to  be  generated  and  saved  ?  1000 
Enter  the  number  of  trials  to  be  generated  ?  5 

Add  a  «gn«i  component  (yM)  [Y1  ?  »  Values  in  square  brackets  [] 

Add  a  clutter  component  (y/n)  ?  are  defatdt  values 

ENTER  PARAMETERS  FOR  THE  CLUTTER. 

Oder  of  the  AR  process  ?  3 

Do  yon  want  tbe  driving  noise  to  be  a  Gaussian  or  SIRP  process  (g/s)  [G]  ? 

Specify  Cbolesky.  L-D*U.  or  SVD  decompositioo  (c/l/s)  [C]  ? 

Do  you  want  tbe  process  to  be  real,  imaginacy.  or  complex  (r^)  [C]  ? 

Use  Gaussian  or  Exponential  shaped  funcdon  (g/e)  [G]  ? 

Enter  tow  1  of  the  amplitude  matrix: 

Emnr  a  row  ?  13) 

Enter  row  1  of  die  one-lag  tenqxiral  oocreladoii  matrix: 

Bner  a  row  ?  0.5 

Enter  row  1  of  the  leg  matrix: 

Emerarow70 

Reference  doppler  frequency  7  0 J 

Enter  the  sample  interval  7  0.01 

A  complex  AR(  3  )  process  will  be  generated. 

The  driving  noise  will  be  Gaussiaa 
Genssiaa  shaped  function. 

Amplitude  matrix:  1.0000e400 
Tenqwral  correladon:  S.OOOOe-01 
Leg  matrix:  0 

Reference  doppler  frequency:  S.OOOOe-01  Hz. 

Sanqtle  intetval:  l.OOOOe-OQ  seconds. 

Are  these  parameters  okay  (yM)  [Y]  7 _ Many  programs  echo  user  inputs. 


Figure  2-2:  Generating  an  AR  Process  with  the  Menu-Based  System 
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GxielidaB  Ug  values: 

R(OH1.0000B400.0.0000e400) 

R(lH4.997Se^)l.lJ7Q5e-02) 

RaH62377e-Q2J.9244e-03) 

R(3M1544Se-03.1.8381»4)4) 

Detanninants  of  priudpie  auuors: 

(1.0000MaO.a000Oe4O0X7.SO0Oe-Ol.0.00Q0e4OOKS:2734e-O1.4.37S7e-lO)(3.65O0e^)1.6.8O6Se-10) 
Coefficients  for  the  AR  process: 
a(l)-  (-6J593e^l.-2.0613e-02) 

aaV  (3^48e-0U.0603e-02) 

a(3)-  (-1  J445e-01.-1.1764e^) 

While  noise  covariance  matrix  (WNCM):  (6.9214e<01.0.0000e400) 

WNCM  Oioies^  decamp.:  (83195e-01.0.0000e+00) 

Add  a  bridle  noiae  component  (yAi)  [Y]  ?  N 

Total  output  file  name  ?  rivar  The  user  saves  the  results  in 

Title  of  this  dataset  ?  An  AR  prooos  a  file  for  later  analysis 

tTuriy  oo^nt  file  ? 

No  file  name  emend,  is  this  okay  ?  y 

Output  coefficient  file  name  ?  rivarxoef 

Tide  of  diis  dataset  ?  AR  Coefficients 

Qntler  AR  coefficients  and  covariances  written  to  file. 

Enter  the  mntdier  of  points  to  be  generated  befom  saving  data  7  5000 
ok? 

MULTI-CHANNEL  (M/C)  PROCESS  SYNTHESIS  MENU 


The  user  returns  to  the  menus 


Figure  2-3:  Generating  an  AR  Process  (Continued) 
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MULTI-CHANNEL  (M/O  DIAGNOSnCS  MENU 

D1  -  DiapUy  Multi-auumel  (M/Q  signal  stats 
D2  -PlatM/Cdau 

D3  -  tafoon  NottaU-Stmid  or  >^dn-Miaf  estiinaiioD 

D4  -  Rafcm  Yule- Walker  esdmatioo 

DS  -- Estunate  All  parameiers  for  each  channel 

D6  -  Esdmaie  covariance 

D7  -  Perfocn  estunatioa  of  stale  space  parameten 

D8  -  Perfcnn  M/C  comlatioa  (teoporal) 

D9  -  Perform  M/C  cooeladoo  (enseinbie) 

DIO  -  Perform  M/C  conelatiao  (quadrature,  temporal) 

Dll  -  Perform  M/C  conelatiao  (quadrature,  ensemble) 

D12  -  Perform  Oztuk’s  algorithm 

D13  -  Plot  eigodic  series  equatiao 

D14  -  (jmenie  periodogram 

DIS  -  Sidit  a  M^  iiqjat  D21  ~  Join  isputs 

D16  -  Add  M/C  signals  D22  -  Subtract  M/C  signals 

D17  -  Convert  data  to  ASQI  file  D23*~  Coovert  ASC3I  file  to  dau 

D18  -  Di^lay  complex  data  D24  -  Display  coefficient  file 

D19  -  Test  modified  Bessel  functiaQ  D25  ~  Calcolate  detectirm  probability 

D20  -  r’aifwiatot  snip  quadratic  form  D26  -  Sum  channels 

LI  -  List  data  files  LZ  ~  show  the  change  log 

*  not  yet  implemented 

Enter  a  commit  or  Q  to  return  to  the  MAIN  menu:  d2 
Tjnking  program  k|dt.  please  stand  by  ... 

'bm^mi4fiqdt’  is  iq)  to  date. 

KFLT  -  omvert  binary  to  text  file  and  plot  under  Openwimlows 
Veisiaal.7 

hqmt  file  name  7  rlvar  The  file  to  plot 

Title:  An  AR  process 

Rest  trial  to  plot  7  1 

Number  of  trials  to  plot  7  1 

X  axis  title  [Hme]  7 

Xaxiaraage  7 

Y  axis  title  (Magnitude]  7 

Data  smoothing  7  0 

Comptex  data  ooitvetsiaa:  Real.  Imaginary.  Magnitude,  or  An^  [R]  7  r 
GOBveiting  data ... 

OiMwwU  written:  1 

Trials  mitten:  1  The  plot  mndow  appears  here 

dt? _ 

Figure  2-4:  Plotting  a  Data  File 
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Magn  i  t  Lide 


An  AR  process 


O  so  100 

T  i  me 


Channe I  1 


Figure  2-6:  The  Final  Plot 
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3.  NEW  CAPABILITIES 

The  ygn^i  detection  algorithm  described  in  Section  2  has  been  implemented  over  a  series 
of  efforts  developing  the  Multichannel  Signal  Processing  Simulation  System.  Typically,  each 
effort  has  implemented  a  new  model  or  a  capability  recently  developed  by  RL/OC  and  their 
contractors.  Under  the  present  effort,  new  capabilities  were  added  to  the  MSPSS  in  the  areas 
of 

•  Sensor  fusion 

a  State  ^ace  models  and  Kalman  filtering 

•  The  Extreme  Value  Method  for  estimating  tail  probabilities 

•  Ozturk’s  algorifiim 

New  capabilities  were  first  implemented  in  the  menu-based  subsystem.  When  appropriate, 
they  were  then  implemented  in  the  UFI-based  subsystem.  This  section  defines  the  capabilities 
implemented  under  this  effmt. 

3.1  Saisor  Fusion 

Signal  processing  algmithms  typically  process  large  volumes  of  data.  These  demands 
often  stress  current  computer  ix'ocessm'  technology.  One  method  for  meeting  this  demand  is 
for  concurrent  processors  to  analyze  data  in  parallel.  A  multichannel  system  could  consist  of 
sensors  for  each  channel  Distributed  processing  would  be  performed  by  each  sensor,  with 
sub^uent  processing  being  performed  at  a  fusitm  center. 

Previously,  the  multichannel  detection  approach  consisted  of  processing  the  observation 
data  vector  using  multichannel  prediction  error  filters  to  obtain  a  whitened  oror  residual. 
These  processes  were  obtained  by  the  simultaneous  removal  of  the  temporal  and  cross-channel 
correlatioa  information.  This  t^roach  required  the  use  of  all  the  multichannel  data  at  the 
processing  center  and  is  referred  to  as  the  centralized  ^>proach.  In  contrast,  the  method 
developed  here  is  a  decentralized  (q)proach  which  individually  {xocesses  data  on  each  channel 
to  first  form  separate  error  signals  rm  each  channel  before  passing  these  signals  to  a  central 
processor  for  whitening  across  chaimels.  The  performance  of  this  scheme  is  suboptimal  to  the 
centralized  apfX’oach,  but  represents  applications  for  which  the  data  is  first  preprocessed  on 
each  channel  (or  sensrH*)  to  reduce  the  redundant  information  content  before  tra^n-al  to  the 
centralized  processor.  In  space-time  signal  processing,  this  case  is  refored  to  as  the  factored 
^)proac^  in  which  temporal  and  spatial  processing  is  performed  distinctly.  Finally,  the 
computational  requirements  of  the  decentralized  approach  are  less  than  the  centralized  case. 

3.1.1  Multichannel  Synthesis  of  AR  Processes 

The  processes  analyzed  by  the  decentralized  approach  are  generated  using  the  existing 
multichannel  MSPSS  synthesis  procedure.  No  new  process  synthesis  routines  were  developed. 
The  radar  return  is  generated  as  the  sum  of  signal,  clutter,  and  noise.  Either  signal  or  clutter 
may  be  absent.  The  noise  is  syntiiesized  as  white  noise,  while  the  signal  and  clutter  are  AR 
processes.  The  sum  of  white  noise  and  AR  processes  is  an  ARMA  process,  which  may  be 
^roximated  as  a  higher  order  AR  process. 

Spedfically,  let  x  ( 1 ),  x  (2),  ...,  x  (/V )  be  a  stochastic  process  ftxmed  by  the  sum  of  three 
processes: 

X  (n)  m  s  (n)  +  c  in)  w(n),  n  ■  1.  2 . N.  (3-1) 
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The  of  vectors  x(n)  represent  the  radar  return,  5(a)  is  the  signal.  c(n)  is  the 

clutter,  and  tv  (a  )  is  the  noise.  Ea^  of  the  subprocesses  can  be  real,  imaginary.  (»-  complex. 
The  number  oi  elements  in  each  vector  x(a),  5(a),  c(a),  and  w(a)  is  the  number  of 
channels  J. 

A  single  sequence  of  vectors  x  ( 1),  x  (2),  ...,  xiN)  is  referred  to  as  a  trial  or  realization 
of  the  stochastic  {u-ocess.  Since  this  is  a  random  process,  it  will  vary  frmn  trial  to  trial.  Each 
vector  X  (n )  in  the  sequence  of  vectors  is  referred  to  as  a  time  sample.  Again,  because  of  the 
random  nature  of  the  process,  the  N  time  samples  vary  across  a  single  realization  of  the 
underlying  process  with  a  q>ecified  temporal  and  o'oss-channel  correlation. 

3  J.1.1  Synthesis  of  Ganasiaii  Noise 

The  noise  vector  w  (a  )  is  synthesized  as  a  zero-mean  Gaussian  white  noise  process  with 
covariance  matrix  ; 

Z^»E[w(n)w"(n)].  (3-2) 

The  noise  process  is  uncocrelated  in  time,  but  may  be  correlated  across  channels  if  the  off- 
diagonal  elements  of  are  nonzero.  By  definition,  !»,  is  Hermitian.  The  components  of 
have  a  physical  interpretation.  The  diagorud  elements  are  the  variances  of  the  corresponding 
channels,  and  the  off-diagoiud  elements  are  the  cross-chaimel  covariances. 

The  user  specifies  either  the  Cholesky  decomposition,  the  LDU  decomposition,  or  the 
Singular  Value  Decomposition  (SVD)  of  the  covariance  matrix.  The  Choleslcy  decomposition 
is  given  by  Equation  3-3: 

(3-3) 

where  C»  is  a  lower  triangular  complex  matrix.  The  LDU  decomposition  is 

(3-4) 

where  is  s  lower  triangular  matrix  with  ones  along  its  diagonal  and  is  s  diagonal 
matrix. 

The  SVD  decomposition  of  I,,  is  given  by  Equatimi  3-5: 

(3-5) 

Aw  is  a  diagonal  matrix  whose  elements  are  ei^values  of  The  columns  of  arc  the 
corresponding  right  hand  eigenvectors,  that  is, 

(3-6) 

The  rows  of  are  the  left  hand  eigenvectors: 

(3-7) 

The  eigenvectors  in  (2*  should  have  a  norm  of  unity.  Also,  since  is  Hermitian,  QH  is  the 
inverse  of  Q,. 

To  goierate  the  white  noise  vectOT,  the  user  specifies  eidier  Lw  and  D^,  or  and 
Aw.  For  each  time  sample  in  a  realization  of  Gaussian  white  noise,  a  vector,  v„(n),  is 
generated  of  statistically  indqwndent  zero-mean  normally  distributed  random  variates.  The  jth 
channel  has  a  variance  of  either  unity,  if  the  user  specified  a  Cholesky  decomposition;  iD„  )j  j, 
if  the  user  specified  an  LDU  decomposition;  or  (Aw  );.y,  if  the  user  specified  a  SVD.  (If  the 
noise  is  cmnplex,  the  real  and  imaginary  components  of  Vw  (a  )  are  indq)endently  generated 
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(3-8) 


with  variances  of  — ,  wbtfe  is  the  desired  variance  of  (n  ).) 

The  white  noise  is  dien  generated  as 

w(n)«r»v»(ii). 

where  is  either  C». or  Q^. 


3J.L2  Synthesis  of  Sigaal  and  Clatter  as  AR  Processes 

The  signal  and  noise  are  each  synthesized  as  Autoregressive  proce»es.  This  section 
describes  the  synthesis  of  the  clutter.  The  signal  is  synthesized  in  an  klenticai  nunner.  The 
AR  process  for  the  clutter  c  ( 1 ),  c  (2) c(N)  is  defined  by  Equaticm  3-9; 

c(«)--  £  A"(it)c(fl -*)  +  £,(n).  (3-9) 

k  m  I 

whoe  c  (n )  is  a  /  X  1  vector  process  and  A/'d),  A" (2) . a"(P  )  are  the  /  x  /  AR  nunix 

coefficients  for  an  AR  process  of  order  P .  £«  (a  )  is  the  driving  noise  term.  It  is  a  zero  mean 
process,  uncorrelated  across  time,  with  covaiiaxwe  matrix  defined  by  Equadtm  3-10: 

I, -£[c(A)c«(n)l.  (3-10) 

The  driving  noise  is  syndiesized  as  white  noise  following  die  procedure  described  in  Section 
3.1.1.!  above.  The  user  ^Mcifies  either  the  Cholesky.  LDU,  or  Singular  Value  Decomposition 
of  the  covariance  matrix  £«.  The  innovations,  unconelated  both  in  time  and  across  channels, 
are  generated  with  ^>propriate  channel  variances.  A  matrix  from  the  approfuiate 
decomposition  is  used  to  premult4)ly  each  time  sample  frmn  the  innovations.  This  process 
yields  the  driving  noise  time  samples  a,  (n )  witii  the  l^)Ixropfiate  conelatitms  across  channels. 

Direct  qtedfication  of  the  AR  coefficients  A^(k )  and  the  covariance  matrix  does  not 
easily  permit  explidt  control  of  the  temporal  and  cross-channel  conelation  of  the  resulting  AR 
process.  Consequently,  Dr.  Michels  has  developed  a  "shaping  function”  approach  for 
contrdling  AR  parameters  in  tenns  of  physically  meaningful  parameters  (Michels  91). 

For  the  stationary  zero-mean  stochastic  process  c  ( 1),  c  (2) . c(N),  a  matrix  showing 

the  correlations  across  channels  for  the  kth  lag  is  defined  by  Equation  3-11: 

R,(l:)-£lc(n)c"(A -*)].  (3-11) 

Note  diat  if  the  number  of  channels  is  /,  R,  (^)  is  a  /  x  /  matrix.  Furthermore,  it  follows 
from  this  definition  that  Equation  3-12  holds: 

£,(-*)-£"(*).  (3-12) 


The  correlation  matrix  and  the  AR  parameters  are  related  by  the  Yule-Walko-  equatimis. 
For  example,  if  there  are  three  lags,  P  =  3,  the  Yule-Walker  equati(»s  are  given  by  3-13; 


[/  A^il)  A^{2)  A" 


R,(0)  R,(l)  RA2)  Pc  (3) 

Pc(-l)  Pc(0)  Pc(l)  Pr(2) 

Pc  (-2)  Pc(-l)  Pc(0)  Pc(l) 
Pc (-3)  Pc (-2)  Pc(-l)  Pc(0) 


Zc  0  0  0 


(3-13) 


Given  tiie  correlation  matrix  Pc.  these  equations  can  be  solved  for  the  AR  coefficients  and  the 
driving  noise  covariance  matrix.  The  NQPSS  uses  the  Levinsmi-Wiggins-Robinson  algorithm 
to  scfive  the  Yule-Walker  equations. 
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The  correlation  matrix  R  is  specified  by  the  user  through  the  use  of  parameters 
describing  certain  shaping  fvmctions  from  which  the  correlation  matrix  is  calculated.  Shaping 
functions  of  a  Gaussian  and  an  Exponential  form  are  available.  The  Gaussian  shaping 
function  for  the  cross-correlation  function  is  given  by  Equations  3-14  and  3- IS  (Michels  90a): 


1  j- 


K.ij 


(3-14) 


(3-15) 


ATe  is  a  cmstant  proportional  to  the  variances  on  each  channel,  A.  is  the  one-lag  temporal 
correlatioa  parameter.  is  the  lag  value  at  which  the  function  peaks.  is  the  Doppler  shift, 
and  r,  is  the  sampling  period,  (x*  is  the  complex  conjugate  of  the  complex  number  x.) 


The  Esqxmential  shaping  function  is  defined  by  Equaticms  3-16  and  3-17; 


1 

^.1.7 


(3-16) 


(3-17) 


3.1.2  The  Multichannel  AR  Model 

The  synthesis  of  signal,  clutter,  and  noise  has  been  described  in  Section  3.1.1.  The  sum 
of  these  processes  is,  in  general,  an  ARM  .  process.  However,  an  ARMA  process  can  be 
approxiinated  by  a  high  order  AR  model.  This  model  may  permit  correlaticm  between 
channels.  Some  of  the  off-diagonal  elements  or  the  AR  coefficients  may  very  well  be  non¬ 
zero.  A  filter  based  on  this  model  would  be  a  tapped  delay  line  filter  with  matrix  coefficients. 
Such  filters  were  fxwiously  implemented  in  the  K^PSS,  but  they  do  not  support  decentralized 
sensor  fusion  analyses. 

Let  x(l),  x(2),  ...,  x(N)  be  a  discrete-time  complex-valued  process.  x(a),  a  =  1,  2,  ..., 
is  a  /  element  vector  where  /  is  the  numbtf  of  channels.  The  multichannel  AR  process  is 
described  by  Equation  3-18: 

J:(«)--  L  A"(*)x(fl )  +  e(n).  (3-18) 

k  m  1 

where  E(n )  is  a  /  x  1  zero-mean  Gaussian  white  (hiving  noise  vector  i»x)cess  with  covariance 
matrix  £.  Since  the  AR  coefficients  a" u )  L  are  matrices,  the  fx-ocess  exhibits  cross- 
channel  correlation. 

The  driving  noise  is  modeled  by  <me  of  three  decompositions  of  £.  The  Cholesky 
decompositicm  is 

I-CC«.  (3-19) 

where  C  is  a  lower  triangular  matrix.  The  corresponding  model  of  the  driving  noise  in 
Equation  3-18  is  given  by  Equation  3-20: 

€(n)«Cv(n).  (3-20) 

where  v(fl)isayxl  zero-mean  unit  variance  Gaussian  noise  vector  uncorrelated  both  in 

time  and  across  channels. 
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The  LDU  decomposition  of  £  is  given  by  Equation  3-21: 

L^LDL**  (3-21) 

D  is  a  diagonal  matrix,  and  £  is  a  lower  triangular  matrix  with  unity  along  the  principle 
diagonal.  The  corresponding  model  of  the  driving  noise  is  given  by  Equation  3-22: 

£(#i)»£v(n).  (3-22) 

In  this  case,  the  covariance  matrix  of  v(n  )  is  D. 

The  SVD  of  £  is  given  by  Equation  3-23: 

£-(2AQ".  (3-23) 

where  A  is  a  diagonal  matrix  whose  elements  are  eigenvalues  of  £.  The  columns  of  Q  are  the 
corre^HMiding  right  hand  eigenvecttvs,  and  the  rows  of  Q"  are  the  left  hand  eigenvectors. 
The  eigenvectors  in  Q  have  a  norm  of  unity.  Also,  since  £  is  Hermitian,  Q"  is  the  inverse  of 
Q.  For  the  S\T>,  the  model  of  the  driving  noise  is  given  by  Equation  3-24: 

e(n  )  ■  (2  v(/i ).  (3-24) 

where  v(fl )  is  a  zero-mean  Gaussian  noise  process  with  covariance  matrix  A. 

Figure  3-1  shows  the  prediction  error  filter  used  to  produce  white  noise  from  a  signal 
synthesized  by  this  model.  This  filter  implements  the  centralized  approach.  The  tapped  delay 
line  structure  whitens  the  ouqput,  e(n).  across  channels.  T  is  either  C,  L,  or  Q  m  the 
Cholesky,  LDU,  or  Singular  Value  Decomposition  respectively.  Thus,  the  final  output  of  this 
filter,  v(a  ),  is  uncoirelated  both  in  time  and  across  ch^ds. 


Figure  3-1:  Tapped  Delay  Line  Filter  for  the  Centralized  Approach 
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3.1^1  A  Diatribated  Tapped  Delay  Lioe  Filter 

In  this  subsection,  we  now  discuss  the  form  of  the  prediction  error  filter  i^ted  in  the 
decentralized  (or  factored)  signal  processing  approach.  Figure  3-2  shows  this  filter.  We  note 
that  this  filter  lacif*  the  off-diagonal  terms  in  the  matrix  coefficients  of  multichannel  prediction 
error  filter  shown  in  Figure  3-1  above. 


Figure  3-2:  Decentralized  Tapped  Delay  Line  Filter  and  Fusion  Center 


The  filter  is  comprised  of  tapped  delay  line  filters  which  first  remove  temporal  correlation 
in  each  channel  j  : 

'‘i 

ej  (n)mxj(n)+  £  a/ffc )xj{n-k),  jm  1. 2 . J  .  (3-25) 

kml 

The  tallied  delay  line  portion  of  the  filters  can  be  distributed  among  the  sensors  for  each 
A  second  phase  of  the  filter,  implemented  at  the  sensor  fusion  site,  removes 
remaining  cocrelatioa  between  channels: 

v(fl)-r*e(ii).  (3-26) 

where  T  is  eidmr  C.  L,  or  (2  in  the  Qiolesky.  LDU,  or  Singular  Value  Decomposition, 
respectively. 

3.U.2  Estiaatioo  of  AR  Coefficieiits 

In  practice,  the  filter  parameters  I  and  aj(k),  k  s  1,  2 . Pj,  j  ~  1,  2 . J,  will  be 

estimates  obtained  from  the  data.  These  estimates  are  obtained  in  two  stages.  First,  the  AR 
coefficients  aj(k )  are  estimated  from  the  jth  chaimel  data.  Then,  the  resulting  error  signals 
from  the  J  chaimels  are  used  to  estimate  the  covariance  matrix  L. 
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Figure  3-3  shows  the  system  structure  used  in  estimating  the  AR  coefficients.  The  input 
is  many  realizations  of  the  modeled  process.  A  set  of  AR  coefficients,  dj il),  a/(2),  .... 
dJiPj ),  is  estimated  for  each  realization  using  some  suitable  single-channel  algorithm.  The 
currently  allows  the  user  to  choose  between  the  Yule-Walker  and  Burg  algorithms 
implemented  as  described  in  (Marple  87).  The  user  must  specify  the  order,  P,,  for  these 
algorithms;  it  is  not  estimated.  These  estimation  algorithms  also  produce  an  estimate  of  the 
variance,  d/  for  each  channel.  These  variance  estimates  are  discarded  and  not  used  further. 


Figure  3*3:  Estimation  of  AR  Coefficients 

Once  a  set  of  estimates  of  AR  coeffidents  is  obtained  for  each  trial,  tiiey  can  be 
averaged  together  across  trials  at  the  user’s  discretion.  This  average  is  the  desired  estimate  of 
the  AR  coefficients  to  be  used  in  the  distributed  tapped  delay  line  filter  shown  in  Hgure  3-2. 
The  sample  variance  of  these  estimates  can  also  be  calculated  across  channels  to  assess 
algoritiun  performance.  .^}pendix  D  describes  a  numerically  stable  one-pass  algorithm  useful 
in  estimating  means  and  variances. 

3.1,ZJ  Eatimalioa  of  the  Covariaiice  Matrix 

The  covariance  matrix  of  the  resulting  error  signals  on  the  /  channels  is  estimated  frcmi 
die  output  data  from  the  /  channels.  As  can  be  seen  firom  Figure  3-4,  die  covariance  matrix 
estimation  process  relies  on  the  previously  determined  AR  coefficients  d/(2),  ..., 

a/  (Pj ).  Ttom  values  can  either  be  actual  values,  if  known,  or  the  estimates  obtained  by  the 
process  described  above.  The  covariance  matrix  is  estimated  from  the  /  x  1  vector  e(ii ) 
produced  by  the  first  phase  in  the  filter. 
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Figure  3^:  Estimation  of  Covariance  Matrix 

The  covariance  matrix  is  estimated  for  each  trial  of  e(A )  as  follows.  First,  the  channel 
mean  is  found  across  time  samples: 

(3-27) 

"  Hal 

where  N  is  the  number  of  time  samples  (die  mean  should  be  statistically  equal  to  zero).  The 
variances  and  covariances  are  then  estimated  by  Equaticm  3-28: 

t.y --i-  £  [ei(«)-i.lU>(/»)-e,]«.  (3-28) 

H  ml 

The  resulting  estimates  are  averaged  across  trials  to  obtain  the  estimate  of  £  whose 
decoaqHmtioo  is  used  in  the  filter. 

3J  Model  Identification 

The  sum  of  noise  and  clutter,  or  of  noise,  clutter,  and  signal,  can  be  modeled  as  a  state 
space  process.  In  the  model-based  signal  detection  algorithm  described  in  Section  2.2,  the 
correyonding  filter  is  a  Kalman  filter.  This  approach  has  been  developed  in  (Roman  93),  and 
the  state  ^ace  model  and  Kalman  filter  were  implemented  in  the  h^PSS  under  this  effort. 
This  implmnentation  provided  t>ie  following  new  capabilities: 

•  C<»version  of  AR  model  parameters  to  state  space  model  parameters 

•  State  space  process  synthesis  based  on  direct  user  input  of  model  parameters 

•  Estifflation  of  state  space  model  parametm 

•  A  Kalman  filter  that  converts  the  radar  return  process  to  Gaussian  white  noise 
(innovations). 
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•  A  signal  detection  algorithm  based  on  the  state  space  model. 

3J.1  A  State  Space  Moitel 

Under  die  state  space  approach,  the  radar  returns  are  modeled  as  follows.  Let  j:  ( 1 ),  x  (2). 
....  X  (N )  be  a  discrete-time  complex- valued  process  representing  the  radar  returns  generated  as 
described  in  Equations  3-29  and  3-30: 

3>(/i+l)»fy(«)  +  Cu(n)  (3-29) 

xin)mH" y(n)  +  D"  w(n).  (3-30) 

X  (rt ),  u  in),  and  win)  are  7 -element  vector,  where  7  is  the  number  of  channels.  The  process 
u  ( 1) . u  iN )  denotes  the  input,  while  w  ( l), ...,  h>  (jv )  is  measurement  noise. 

yin)  denotes  an  Af -element  vector  rqiresenting  the  state  variables.  The  state  space  can 
be  diought  of  as  a  mem(X7  associated  with  the  process.  The  state  variables  evolve  in 
accordance  with  the  first  equation  above,  while  the  value  of  the  output  process  at  any  time  is 
a  transformation  of  the  state  variables  and  measurement  noise. 

3,2.L1  The  Innovatioas  Representation  of  the  State  Space  Model 

To  implement  an  innovations-based  signal  detection  method,  filters  are  constructed  such 
that  the  outputs  of  the  filters  are  the  innovations  under  the  tqipropriate  hypotheses. 
Innovations  are  defined  by  a  (not  necessarily  Gaussian)  white-noise  process.  Kalman  filters  are 
used  to  produce  the  innovations. 

The  state-space  model  can  be  transformed  into  the  following  innovations  representation: 

a(« +  l)«Fo(n)  +  A:6(/i)  (3-31) 

xin)mH^  ain)-¥ein).  (3-32) 

The  process  ed),  e(2),  ...,  ziN)  defines  the  innovations,  and  the  process  a(l),  ail),  ...,  aiN) 
is  the  state  of  die  innovations  rq)resentation.  The  matrix  AT  is  the  Kalman  filter  gain. 

Figure  3-5  shows  a  system  block  diagram  for  the  state  space  model.  Although  the  output 
process  and  certain  matrix  parameters  are  r^resented  by  the  same  symbols  as  in  Equation  3- 
29  and  3-30,  these  parameters  may  reflect  a  change  in  basis,  liie  innovations  model  is 
"correlatitm  equivalent"  to  the  general  state  space  model. 

32.12  A  Basis  Transfonnation  of  die  laaovatioes  Represeolation 

The  innovations  rq)resentation  of  the  state  space  model  is  defined  tmly  up  to  a  basis 
transformatitm  of  the  state  vector.  Let  r  be  an  invertible  M  x  M  matrix  representing  a  basis 
transformation  of  the  state  vecttx-.  Define  fl(f( ),  Kp,  and  as  follows: 


P(n)*r ain ) 

(3-33) 

Fp-rFT-' 

(3-34) 

KpmTK 

(3-35) 

(3-36) 
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Figure  3-5:  Innovations  Representation  of  State  Space  Process  Synthesis 

Under  this  basis  transformaticm,  an  alternative  innovations  rq)resentation  of  the  state  ^ace 
model  is  pven  by  Equadcm  3-37  and  3-38: 

p(«  (3-37) 

*(«)-«{fp(ii)  +  e(H).  (3-38) 

The  innovations  and  the  output  process  remain  unchanged. 

Equations  3-37  and  3-38  are  of  the  same  form  as  Equations  3-31  and  3-32.  Both  sets  of 
equations  are  innovations  rq)resentations  of  the  same  state  space  model.  This  demtmstrates 
-  tfam  the  innovations  rqnesentation  is  unique  only  up  to  a  basis  transformation  of  the  state 
space.  A  practical  implication  is  a  difficulty  in  validating  state  ^Moe  estimation  algorithms. 
The  output  process  may  be  synthesized  with  given  state  qMce  parameters,  but  estimates  of 
diese  parameters  obtained  from  the  omput  process  can  dUfo*  from  tfae  parameters  used  in 
syndi^.  These  differences  sre  acoeptiU)le  as  long  as  they  represent  a  baais  transformation  of 
die  original  parameters.  Whether  this  is  so  or  not  may  be  difficult  to  determine.  In  any  case, 
the  model  oi^  and  the  covariance  matrix  for  the  innovations  sre  unique. 

32.L3  A  Slate  Space  Represeatetioa  of  Aa  Aatoretreanve  Model 

An  Autoregressive  (AR)  model  can  be  cast  in  the  form  of  a  state-space  model.  Consider 
the  case  of  the  /*th  order  AR  model  x  ( l ),  x  (2), ...,  x  (N ): 

x(rt)-- £  A" (t)x(« -*)  +  £(«).  (3-39) 

*  >1 

where  eCn )  is  a  zero  mean  driving  noise  term  with  covariance  matrix  £. 
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In  an  AR  model,  the  current  value  of  the  process  is  a  linear  combination  of  the  last  P 
values  and  white  noise.  In  a  state  space  model,  the  current  value  is  a  function  of  the  state 
variables  and  white  noise.  Thus,  let  the  state  variables  be  the  past  P  values  of  the  process: 


'  x(n-l) 
x(n  -2) 


Lx(« -P)i 


(3-40) 


Since  x  (a  )  is  a  /-element  vector,  there  are  A/  =  /*  x  /  state  variables. 


To  ccmiplete  the  state  space  representation  of  the  AR  process,  let  the  state-space 
parameters  be  defined  as  follows: 


-A"(l) 

-A"  (2) 

-  1) 

-A«(P) 

/ 

0 

0 

0 

0 

/ 

0 

0 

0 

0 

/ 

0 

(3-41) 


(3-42) 


H"  m[-A"  (2)  •  -A"(P)] 


(3-43) 


D»ml  (3-44) 

Furfiiermore,  let  the  ii^)ut  noise  and  the  measurement  noise  be  equal  to  the  driving  noise 
{M-ocess  in  the  AR  model: 

w(H)mu(n)mt(n).  (3-45) 

Using  fiiese  definitions  of  the  relevant  parraetm.  the  iiq>ut  process,  and  the  measurement 
ndse,  die  AR  process  is  rei»esented  in  a  state  space  form,  as  defined  in  Equations  3-29  and 
3-30.  This  state  space  rq>resentation  is  also  the  iimovations  representation  of  the  state  space 
model. 


3J.L4  A  State  Space  Represeatatioa  of  Aa  ARMA  Model 

An  Autoreghessive  Moving  Average  (ARMA)  model  can  also  be  conceptualized  as  a 
state-space  model.  Define  the  ARMA  model  to  be  generated  by  the  sum  of  AR  and  white 
noise  (vocesses: 

x(n)ms(n)*ei(n).  (3-46) 

where 
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i  («)--  £  M"(i)5(n -*)  +  e2(n).  (3-47) 

i  m  I 

This  is  s  special  type  of  ARMA  model,  and  the  state  space  rq)resentation  defined  below  for 
tfiia  model  is  not  valid  for  all  types  of  ARMA  models. 

In  the  state  space  rcfnesentadon,  the  state  variables  are  the  M  =  P  x  J  variables  defined 
by  past  values  of  the  AR  process; 

i(ii  - 1) 

s(n-2) 

(3-48) 

s(n-P). 

The  state  space  parameters  are  defined  as  in  the  AR  model; 

-A^il)-A^(2)  ■■■  -4"(/>  -  1)  -4''(/») 

/  0  0  0 
0/0  0 
f  -  .  (3-49) 


0  0/  0  J 

/ 

0 

(3-50) 

0 

//" -[-4"(l)-4«(2)  -4''(P)]  (3-51) 

/)"-/  (3-52) 

The  state  space  representation  of  this  type  of  ARMA  process  is  distinguished  firom  an 
AR  process  by  the  treatment  of  the  iiq>ut  process  and  die  measurement  noise; 

u(n)mt2(n)  (3-53) 

»*'(«)-ei(ii)  +  e2(n)  (3-54) 

Notice  that  the  resulting  state  space  model  is  not  an  innovations  representation.  Conversion  of 
this  model  to  an  innovations  representation  requires  determinatitm  of  the  Kahnan  gain  fw  this 
system. 

Other  special  cases  of  the  state  space  model  are  also  of  interest,  although  not  defined 
here.  For  example,  the  multipath  model  defined  in  (hficbels  91)  can  be  rq>resented  by  a  state 
space  model. 
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Innovatiofis  Generation 

Given  a  state  space  model  of  the  process  xin),  a  Kalman  filter  can  be  used  to  generate 
the  innovations  sequence.  Alternatively,  the  whitening  filter  associated  with  the  Kalman  filter 
can  be  used.  The  MSPSS  uses  a  ste^y  stam  approximation  to  a  one-step  predictor  Kalman 
filter  followed  by  a  linear  transformation.  This  filter,  is  shown  in  Figure  3-6. 


Figure  3-6:  The  Innovations  Generation  Filter 

Formally,  let  x  ( l ),  x  (2) . jc  (/v )  be  the  input  process  to  the  filter.  The  first  phase  is  a 

Kalman  filter,  and  its  output  is  the  process  e(l),  e(2),  ....  t{N ).  The  second  phase  is  a  spatial 
filter  which  removes  ctxrelation  between  chaxmels.  The  output  of  this  second  phase,  and 
therefore  the  ouq}ut  of  the  entire  filter,  is  v(  l ),  v(2} . viN ). 

The  ouq)ut  of  the  first  phase  at  each  time  sample  is  the  difference  between  the  input  and 
a  minimum  variance  estimate  of  the  input; 

e(«)»ac(n)-i(nln-l).  (3-55) 

The  estimate  of  the  input  is  based  on  an  estimate  of  the  state  variables  in  the  iimovations 
rqvesentation  of  the  state  space  model; 

x(n  In  -  1)»A/"  d(n  In  -  1).  (3-56) 

The  estimate  of  the  state  variables  is  updated  in  accordance  with  the  system  dynamics: 

d(n  +  1 1  n  )  ■  F  d(n  I  n  -  1)  + /f  e(n  )  (3-57) 

d(0l-l)-0,  (3-58) 

where  /T  is  the  "Kalman  gain  matrix,"  and  the  matrix  F  embodies  the  "system  dynamics."  The 
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product  of  the  gain  matrix  and  the  output,  Kein),  is  known  as  the  "correction  vector."  In 
many  filter  implementations,  the  gain  matrix  is  updated  for  each  time  sample.  The 

fact  that  1C  is  not  updat^  for  this  Kalman  filto*  reflects  a  steady  state  ^}proximati(m. 

The  second  phase  of  the  filter  is  based  on  a  decomposition  of  the  covariance  matrix  L  of 
the  driving  noise  process  in  the  innovations  rqxesentation  of  the  state  space  model.  Let  r  be 
either  C,  L,  or  Q  in  the  Cholesky,  LDU,  or  Singular  Value  Decomposition  of  L,  respectively. 
The  ou^ut  of  die  second  phase  is  calculated  by  decorrelating  the  ou^ut  of  the  first  phase 
across  channels: 

v(fl)-r-‘£(fl).  (3-59) 

When  a  SVD  of  £  is  used,  r~‘  is  equal  to  T" .  Hence,  the  inverse  of  Q  is  calculated  as  its 
Hermetian  transpose. 

If  the  filter  parameters  embody  the  state  space  model  firom  which  the  input  is  generated, 
the  ouqiut  of  this  Kalman  filter  will  be  a  reasonable  estimate  of  the  innovations  process.  In 
practice,  the  true  filter  parameters  are  unavailable,  and  estimates  are  used  for  the  matrices  K, 
F,  H",  and£. 

3,2.3  Estimation  of  State  Space  Parameters 

The  MSPSS  can  be  used  to  synthesize  AR  and  ARMA  processes,  as  well  as  a  few  more 
complex  processes  (e.g.  for  multipath).  As  shown  in  Section  3.2.1,  AR  and  ARMA  processes 
can  be  rqiresented  exactly  in  a  state  space  model.  This  section  defines  the  Canonical 
Correlations  Algorithm,  which  is  used  to  estimate  state  space  model  parameters  in  the 
MSPSS. 

Let  X  ( 1),  X  (2),  ...,  X  (A/ )  be  a  single  realization  of  a  stochastic  process  generated  by  a 
state  q)ace  model.  The  Canonical  Cmrelations  Algorithm  uses  data  on  Ng  such  realizations  to 
generate  an  estimate  of  the  state  space  model  order  M  and  the  matrices  f ,  K,  and  A/"  which 
tq^^ear  in  the  innovations  rejnesentatitm  of  the  state  space  model  (Section  3.2.1. 1).  The 
covariance,  £,  of  the  innovatitms  is  also  computed. 

3,23.1  Step  1:  Eatiwate  Corrdatioa  Matrices 

The  algoritiun  is  based  on  estimates  of  correlation  matrices,  not  the  raw  data.  These 
matrices  are  easily  calculated.  Let  A/  represent  the  correlatimis  both  within  a  channel  and 
between  channels  at  lag  /: 

A, -£[x(«)x"(n -/)].  (3-60) 

Since  x(n)  is  a  /-element  column  vector  and  its  Hermitian  transpose  is  a  /-element  row 
vector.  A/  is  an  array  with  /  rows  and  /  columns.  A,  is  a  Hennitian  matrix  only  for  the  zeroth 
lag.  The  Hermitian  transpose  of  correlation  matrices  for  positive  lags  provide  a  convenient 
metiiod  for  calculating  correlations  for  negative  lags: 

Aw -A,".  (3-61) 

Given  data  on  Ng  realizations  of  a  stochastic  process,  the  correlation  matrices  can  be 
estimated  by  averaging  the  sample  correlations  over  time  and  across  all  trials.  Two  methods 
exist,  one  producing  biased  estimates  of  the  correlation  matrices  and  the  other  producing 
unbiased  estimates.  For  nonnegative  /,  the  biased  estimate  is  generated  as 
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I  x'(l:){x'(4:-/)f  .  (3-62) 

.  .1  "  ik  .{■*.! 

where  x'(n)  is  the  nth  time  sample  of  the  vector  x(n)  from  the  ith  realization  of  the 
stochastic  process.  The  imbiased  estimate  is  generated  as 

i  x‘(*)(x‘(it-/)J«.  (3-63) 

-  ‘  *-i>i 

The  correlation  matrices  Ao,  Ai . At  _  i  are  used  to  estimate  the  past  and  future  block 

correlatitm  matrices.  The  estimate  of  the  /  L  xJL  past  correlation  matrix,  in  block  form,  is 


Rp-.L.L^ 

Ao  Ai 

A»i  Ao 

k-i 

■  k-2 

x(/i-l) 
x(fl  -2) 

lx"(n-l)  x»(n-2)  .  .  .  x"(«-L)l  (^-64) 

Ai-i,  Aj.i, 

Ao 

xin-L). 

The  estimate  of  the  /  L  x  7  L  future  block  cmrelation  matrix  is 


Ao  A_i 

A]  Ao 

k-L 

^1-L 

•E 

x(n) 
x(n  + 1) 

Ix^'in)  x"(«+l)  .  .  .  x"(n+L-l)] 

1 

1 

Ao 

Mfi-t-L-l). 

Notice  that  die  future  block  conrelatitHi  matrix,  Rf-.l.l*  is  the  block  transpose  of  the  past  block 
correladcm  matrix,  Furthermore,  both  the  past  and  future  block  correlation  matrices 

are  Hdmitian. 

The  ccmstant  L,  which  defines  the  number  of  block  rows  and  columns  in  the  past  and 
future  block  correlation  matrices,  is  a  user  input.  It  should  be  between  one  and  20,  and  chosen 
such  that  the  user  believes  /  L  to  be  greater  than  M,  the  number  of  state  variables. 

32J2  Step  2:  Cakulate  Square  Root  Matrices 

Square  root  matrices  are  calculated  for  the  past  and  future  block  conelation  matices.  The 
square  root  is  obtained  using  the  Singular  Value  Decomposition  (SVD)  of  the  block 
correlation  matrices.  The  SVDs  of  the  past  and  friture  block  conrelation  matrices  are  defined 
by  Equations  3-66  and  3-67: 

RF..L.fUF5FUp.  (3-66) 

Rf:l.l^VfSfVP.  (3-67) 

The  matrices  Sp  and  Sf  are  diagonal  nuitrices  composed  of  the  eigenvalues  of  the  block 
correlatitm  matrices  arran^d  in  decreasing  order  of  magnitude.  The  eigenvalues  of  the 
matrices  in  Equations  3-64  and  3-65  are  positive.  The  square  roots  of  matrices  Sp  and  Sp  can 
be  calculated  cm  an  element-by-element  basis  along  the  princq)al  diagonal.  Since  the  block 
correlation  matrices  are  Hermitian,  the  inverses  of  Vp  and  Up  are  their  Heimitian  transposes. 
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Therefore.  ±e  following  hold: 


Rf.L.L^iUfS^UifnUrS^Uif).  (3-68) 

Rf -.lx  - ( Uf UP).  (3-69) 

Thus,  the  square  roots  of  the  block  correlation  matrices  can  be  calculated  as  shown  in 
Equations  3-70  and  3-71: 

^F?L.fU^Sf^Up.  (3-70) 

Rf'^lx^UfSP^UP.  (3-71) 

3JL3J  Step  3:  Cakolate  The  Intermcdtate  Matrix  A 

An  estimate  of  the  /  L  xJL  stochastic  block  Hankel  matrix  is  found  as  follows: 


Ai  Aj 

Aj  As 

■  aJ 

- 

4x,4.i 

x(n) 

x(n+l) 

1 

Al  At+i 

■  Ml  _  1 

mE 

.x(/i+L-l). 

[jc"(/i-l)  x"(rt-2)  .  .  .  x^(n-L)J  0-’^’^) 

The  intermediate  J  L  xJ  L  matrix  A  is  defined  as  follows: 

^  -  Rf'fi.L  Hl.l  Rf'ftx  -  ( Vf  Sf^ UP)Htx  ( U,  S,-'« UP).  (3-73) 

3  JJ.4  Step  4:  Parftw  the  SVD  of  A 

An  SVD  of  A  is  performed  next.  The  SVD  is  defined  as: 

A^U^S^VP,  (3-74) 

wtoe  Sa  is  a  diagtmal  matrix  with  the  singular  values  of  A  along  the  principal  diagonal  in 
order  of  decreasing  magnitude.  All  the  singular  values  are  positive  r^  numbers  between 

.  zero  and  unity.  Let  p,.  P} . pyx,  denote  the  (oxlered)  singular  values.  These  singular  values 

are  the  correlations  between  the  normalized  past  and  the  normalized  future,  and  are  referred  to 
as  foe  canonical  correlafions. 

3.2JL5  Step  5:  KstiaMte  Model  Order 

The  model  order  (foe  number  of  state  variables),  M,  can  be  determined  firtmi  the 
canonical  correlations.  In  theory,  exactly  M  cantmical  correlations  should  be  nonzero,  but  due 
to  numerical  and  statistical  estimation  errors  some  may  be  nonzero.  Methods  fcx  determining 
model  order  are  based  on  identifying  the  number  of  canonical  correlations  that  are 
meaningfully  diffoent  from  zero. 

A  number  of  methods  exist  for  determining  the  model  order.  These  methods  are  based 
on: 

•  A.  The  canonical  correlations 

•  B.  The  normalized  nmning  sum  of  the  canonical  correlations 
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•  C.  The  sqiiare  of  the  cancmical  correlations 

•  D.  Tire  normalized  running  sum  of  the  squares  of  the  canonical  correlations 

•  E.  Log  parameters 

•  F.  Normalized  mutual  informatitm  parameters 

lire  user  can  drcument  the  whole  {M-ocess  of  calculating  the  model  order  by  entering  the  order 
direcdy. 

Affethod  A. 

The  model  order  M  can  be  determined  as  the  largest  value  such  that  for  all  >  M.  the 
canonical  correlation  pt  is  less  than  or  equal  to  scant  user-specified  (small)  value  between 
zero  and  unity. 

Method  B. 

The  normalized  running  sum  of  the  canonical  correlations  is  defined  to  be; 

i'pj 

NRS(k)m  ^ - .  (3-75) 

£  lp.  l 

<  >  1 

In  theory,  all  values  of  p,  should  be  nonnegative,  but  numerical  errors  may  lead  to  small 
negative  values.  The  model  order  can  be  selected  as  the  smallest  value  M  such  that  for  all  k 
2  A#,  NRS  (k)  is  greater  than  or  equal  to  some  user-specified  (large)  value  between  zero  and 
unity. 

Method  C. 

The  squares  of  the  cantmical  correlations  are  merely  p^  pi . p/l-  The  model  order  M 

can  be  determined  as  the  largest  value  such  diat  for  all  Jt  >  Af .  the  canonical  correlation  Pt  is 
less  than  or  equal  to  some  user-qredfied  (small)  value  between  zero  and  unity. 

Method  D. 

The  normalized  running  sum  of  the  squares  of  the  canonical  cmrelations  is  defined  to  be; 

k 

NRSS(k)m^j^.  (3-76) 

Zp? 

i  m  I 

The  model  order  can  be  the  smallest  value  M  such  that  for  all  it  2  Af.  NRSS  ik )  is  greater  than 
or  equal  to  some  user-specified  (large)  value  between  zero  and  unity. 

Method  E. 

The  log  parameters  are  defined  as  -in(l-p|}.  -ln(i-p|).  ....  -ind-p/^).  Since  the 
cantmical  correlations  are  sorted  in  decreasing  order,  the  log  parameters  are  also  in  decreasing 
order.  The  model  order  M  can  be  determined  as  the  largest  value  such  that  for  all  it  >  A/,  the 
log  parameter  -]n(l  -  p|)  is  less  than  or  equal  to  smne  user-specified  (small)  value  between 
zero  and  unity. 
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Metliod  F. 

Tbe  noraulized  mutual  informatioa  parameters  are  defined  to  be: 

i  tad -Pi*) 

-  (3-77) 

2:  ta(l-p?) 

i  >t 

The  model  order  can  be  the  smallest  value  Ml  such  that  for  all  Jk  2  Af.  is  greater  than  or 
equal  to  some  user-specifi'sd  (large)  value  between  zero  and  unity. 

Model  orcter  determinatitm  is  an  important  step  in  tbe  canonical  correlatitms  algmithm, 
and  can  be  implemented  using  any  of  the  methods  defined  above. 

32J.6  Step  6:  Coapate  TraasCormalkia  Matrices 

The  two  JL  X  JL  "transformation''  matrices  Tp  and  Tf  are  defined  by  Equations  3-78 
and  3-79: 

Tp  -  (3-78) 

TfU^Rf'i.L  (3-79) 


3^.7  Step  7:  Eadaute  the  System  Matrix 

The  column-shifted  JL  xJ  L  block  Hankel  matrix  is  defined  as  follows: 


Aj 

A3 

■  ■  ■  Ai.,1 

A3 

A4 

•  •  •  At.,2 

At*  I 

^  +2 

•  Aji 

Next  tbe  following  JL  xJL  matrix  is  calculated: 


(3-80) 


(3-81) 

Let  Z  be  die  square  matrix  formed  fi’om  A  by  taking  the  first  M  rows  and  columns,  where  M 
is  the  model  order  detennined  using  one  of  tbe  methods  sped&ed  in  Stq>  5.  Aside  from 
numerical  errors,  all  columns  of  A  beyond  dm  ^h  column  are  zero. 


Tbe  system  matrix  F  is  estimated  as: 

F-SA.fZSx.f. 


(3-82) 


where  is  the  M  x  3/  square  matrix  in  die  upper  left  of  obtained  using  the  SVD  of  A 
as  defined  in  Stq>  4.  Matrix  F  is  the  system  matrix  in  the  innovatitms  rqiresentation. 


32.3JS  Step  8:  Estimate  Obscrvatioa  Matrix 

Let  Zii  denote  the  matrix  formed  firmn  the  first  J  rows  and  M  columns  of  the  matrix 
Hix  fP-  'ITien  the  Hermitian  transpose  of  the  observation  matrix  H  is  estimated  as  follows: 
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(3-83) 


This  is  the  observation  matrix  in  the  innovatioi^  representation. 

3 JJ.9  Step  9:  Estimate  Backward  Modd  Observatioo  Matrix 

Let  Zr  denote  the  first  M  rows  and  J  columns  of  the  matrix  Tf  ■  ’^cn  die  backward 
M  xJ  observation  matrix  r  is  estimated  as  follows: 

r-SA.'PZr.  (3-84) 

3.2J.10  Step  Mh  Determiae  the  Remaiaing  System  Model  Parameters 

The  remaining  system  parameters  are  estimated  as  follows: 

n-Sx.i  (3-85) 

(3-86) 

K»[T-FnHVr\  (3-87) 

In  the  innovations  representation,  n  denotes  the  (zeroth  lag)  correlation  matrix  of  the  state 
space  vector,  d  denotes  the  (zeroth  lag)  conelation  matrix  of  the  multichannel  innovations 
process,  Ao  denotes  the  (zeroth  lag)  ctxrelation  matrix  of  the  model  output  process,  and  K 
denotes  die  Kalman  gain. 

3J  Extreme  Value  Method 

The  determination  of  thresholds  for  given  false  alarm  probabilities  is  one  amiability  of 
die  MSPSS.  Thresholds  are  found  by  simulating  many  realizations  of  the  stochastic  process 
corresponding  to  the  null  hypothesis,  which  omtains  no  signal.  These  realizations  are  passed 
dirou^  the  signal  detection  algorithm,  thovby  yielding  a  list  of  estimates  of  the  loglikelihood 
statistic.  A  threshold  is  estimated  as  a  value  of  the  statistic  such  that  the  probability  of 
exceeding  this  threshold  is  the  given  false  alann  jx-obability.  Since  the  false  alam  probability 
is  a  tail  {xobability,  an  accurate  estimation  oS  thresholds  requires  the  generation  of  a  very 
large  number  of  realizations  of  the  modeled  stochastic  process.  A  full  simulation  analysis  of 
a  signal  detection  algorithm  is  thus  quite  expensive  in  terms  of  time. 

Prakosh  Chakravarthi  (92)  has  developed  a  method,  known  as  the  Extreme  Value  Theory 
(EVT),  to  accurately  estimate  thresholds  based  <m  an  order  of  magnitude  less  realizations,  and 
this  method  was  implemented  in  the  MSPSS  under  this  effort.  The  EVT  im>proximates  the  tail 
of  the  distribution  of  the  loglikelihood  statistic  by  a  Generalized  Pareto  distribution.  The 
parameters  of  the  Generalized  Pareto  distributicm  can  be  estimated  with  fewer  samples  of  the 
lo^ikelihood  statistic  than  is  needed  to  directly  estimate  a  threshold.  Once  the  C^neralized 
Pareto  distribution  is  fiilly  known,  thresholds  can  be  calculated  for  any  given  false  alarm 
probabilities.  The  programs  for  estimating  thresholds  were  modified  to  permit  the  user  to 
choose  this  method. 

Let  ^  A^^’  ^  ...  ^  A^*’  be  ordaed  sample  statistics  for  the  loglikelihood  ratio  generated 
under  the  null  hypothesis  without  a  signal  present.  Let  F  be  the  underlying  cumulative 
probability  distribution  fw  the  (unordered)  loglikelihood  statistic,  and  let  p  draote  the  desired 
false  alarm  probability.  The  problem  treated  here  is  to  determine  an  estimate,  il,  of  the 
threshdd  such  that 
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(3-88) 


where  is  an  estimate  of  n- 

The  EVT  (Rangaswamy  93)  is  based  on  aj^roximating  the  tail  of  F.  where  the  tail  is 
defined  to  be  the  following  set  of  loglikelihood  ratios: 

{XiXiAo.  l-F(>«)-a}.  (3-89) 


where  a  is  a  user-defined  (smaU)  probability.  (In  other  words,  Pr(A^1<o)^a.)  The  tail  of 
almost  any  distribution  F  can  be  tqiproximated  by  the  Generalized  Pareto  Distribution.  The 
distributimi  function,  C ,  for  the  Generalized  Pareto  Distribution  is  given  by  Equation  3-90: 


G(X) 


1- 


o 


-Ur 


(3-90) 


The  tail  of  F  is  approximated  as  follows: 

P(A.)-(J -a)  +  aG(3l-Ao).  A.>Xo.  (3-91) 


Given  estimates  of  Ao,  y,  and  o,  an  estimate  of  the  appropriate  threshold  can  be  found  for  any 
false  alarm  probability: 


■  Ao+ 

y 


(3-92) 


(This  formula  disagrees  with  Equation  10.51  in  (Rangaswamy  93),  which  is  mistaken.) 

Once  procedures  for  estimating  the  tail  pot^tile  and  the  parameters  of  the  Generalized 
Pareto  Disttibutimi  are  given,  the  EVT  is  completely  specified.  The  tail  percentile  is  easily 
estmuMd.  Let  /  -[a/v.J,  where  [a/V.J  i,  the  imeger  pm  of  «hr,.  Then,  the  foUowing 
equatitm  can  be  used  as  an  estimate  of  the  tail  percentile: 

(3-93) 

The  tail  is  defined  to  be  those  loglikelihood  ratios  A*'’  such  that  A<'^  >  Ao.  If  is  tied  with 

succeeding  values,  the  tail  may  contain  less  than  j  values.  Let  Z(])  ^  Z(2)  ^  ...  s  Z^„y  denote  the 
tail  values,  that  is 


7  _  A*'''*  a 

Z(,)»A  -Ao- 


(3-94) 


Three  methods  are  available  for  estimating  the  parameters  of  the  Generalized  Pareto 
Distribution: 

•  Maximum  Likelihood 

•  Probability  Weighted  Measures 

•  Ordered  Sample  Least  Squares 


3J.1  Maximum  Likdihood 

The  maximum  likelihood  method  is  based  on  reparameterizing  the  likelihood  function  for 
the  tail  statistics  in  terms  of  a  and  t,  where 

T-T/a.  (3-95) 

The  derivative  of  the  reparameterized  loglikelihood  function  with  respect  to  o  is 
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(3-96) 


. ^  Iii(  1  +  ). 

da  a  (tt  ; 

(Equation  10.22  in  (Rangaswamy  93)  contains  an  error  in  transcripti<Hi.)  A  value  that 
the  loglikelihood  ratio  can  be  found  by  setting  this  derivative  equal  to  zero.  The 
estimate  for  o  is  then  found  to  be: 

a(^)--^  £  lo(l  +  Tz.  ).  (3-97) 

m  T  i  .1 

This  expression  can  be  used  to  express  the  rq)arameterized  loglikelihood  fiinctitm  in  terms  of 
t  alone: 

-lnZ,(T;r, . r,)«mlna(T)+  1  +  ~  5^1n(l+Tz,).  (3-98) 

0(1 )t  i  ■ 1 

An  estimate  for  -c  is  obtained  by  minimizing  the  above  equation.  The  system  minimizes  this 
function  numerically  by  using  the  Nelder-Mead  algorithm.  The  implementation  of  the  Nelder- 
Mead  algorithm  in  the  MSPSS  is  based  on  that  in  (Press  86).  This  numerical  method  requires 
two  initial  guesses  fa*  t.  Our  implementation  uses  two  inital  guesses  based  (m  the  Probability 
Weighted  Moments  estimates,  which  can  be  found  in  closed  form.  One  is  three  halves  of  the 
Probability  Weighted  Moments  t,  and  the  other  is  one  half  the  Probability  Weighted  Moments 
1. 

Once  i  and  d(T)  are  found,  y  is  estimated  by 

Y-6(f)^.  (3-99) 

The  above  {x-ocedure  is  obtained  from  (Rangaswamy  93).  We  did  not  check  that  the 
likelihood  fimcticm  is  indeed  maximized  by  diese  solutions  and  that  the  estimates  so  obtained 
do  not  characterize  a  minimum  or  a  saddle  ptwt. 

3 J J  Probabflitj  Weighted  Moments 

Let  Eo  and  e,  be  the  following  probability  weighted  moments  of  the  Generalized  Pareto 
Distributitm: 

eo-£[Z).  (3-100) 

and 

e,-£(Z(l-C(Z))].  (3-101) 

These  parameters  are  estimated  from  the  data  as  follows: 

‘•-i.l,--'  (3-102) 

Calculating  die  relevant  integrals  yields  equadcms  for  the  probability  weighted  moments 
in  terms  of  the  parameters  of  the  Generalized  Pareto  DistributitMi: 
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a 


(3-104) 


Co 


1-Y 


®‘"  2(2-y) 


(3-105) 


These  equations  can  be  solved  for  the  Generalized  Pareto  Distribution  parameters.  The 
following  estimates  of  the  parameters  are  thereby  obtained: 


2loei 

I0-2I, 


(3-106) 


(3-107) 

lo  Io--2ei 

(Equation  10.32  in  (Rangaswamy  93)  is  another,  but  equivalent  formulation.) 


3  J.3  Ordered  Sample  Least  Squares 

The  derivatitHi  of  the  Ordered  Sample  Least  Squares  method  is  given  in  (Rangaswamy 
93).  The  e:q)ected  value  of  each  order  statistic  in  the  tail  can  be  found  as  fimction  of  y  and  a. 
The  estimators  of  these  parametm  are  chosen  to  minimize  the  sum  of  the  squares  of  the 
differences  between  the  observed  order  statistics  and  their  expected  values. 

Using  this  procedure,  an  analytical  fcxmula  can  be  found  for  o  in  terms  of  y  and  the 
observed  data: 


E  r,(Q,(Y)-l) 

d(Y)-Y^^r - -  (5*108) 

E  (Qr(Y)-l)* 

r-l 

where 

(3-109) 

The  parameter  y  is  found  by  a  numerical  fx-ocedure  that  minimizes  the  following  function: 

5(y)- E  (^r--^[Qr(Y)-ll)*  (3-110) 

We  did  not  check  that  a  solutitm  exists  to  tiiis  minimization  problem.  The  numerical  method 
requires  two  initial  guesses  for  y-  Our  inq)lementatioo  uses  two  inital  guesses  from  the 
Probability  Weighted  Moments  estimate  of  y-  These  guesses  are  three  halves  and  oat  half  the 
Probability  Weighted  Mtxnents  y. 

3.4  Ozturk*s  Algorithm 

Dr.  Aydin  Ozturk  has  invented  a  powerfril  algorithm  for  determining  the  distribution  of  a 
random  sanq>ie  (Ozturk  90a,  90b,  and  90c).  It  provides  a  graphical  rqxesentation  of  data.  The 
algorithm  can  be  ^}plied  in  two  modes.  Under  the  first  mode,  it  {vovides  a  formal  statistical 
hypothesis  test  of  the  goodness  of  fit  of  a  sample.  The  other  mode  provitfes  a  means  of 
estimating  the  distribution  of  the  sample.  These  modes  contrast  with  other  goodness-of-fit 
tests,  which  usually  do  not  provide  an  indication  of  an  appropriate  distribution  if  they  reject 
the  null  hypothesis,  where  the  null  hypothesis  is  that  the  data  come  from  a  specified 
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disthbutiim.  Both  modes  were  implemented  under  this  effmt. 

Typically,  Ozturk’s  algorithm  cannot  be  applied  to  the  data  synthesized  by  the  MSPPS  in 
its  raw  The  MSPSS  synthesizes  multichannel  complex  data,  while  Ozturk’s  algorithm 
applies  to  real  univariate  data.  Hence  the  probto  arises  concerning  how  to  apply 
Ozturk’s  axgcrithm  to  cmnplex  data.  Two  possibilities  exist.  Ozturk’s  aigoridim  can  be 
extended  to  multichannel  data  (Romeu  90).  Or  procedures  could  be  created  to  convert 
multichannel  cmnplex  data  to  univariate  real  data.  The  latter  approach  was  adopted  fw  this 

t^if 

3.4.1  Splitting  Channels 

The  simplest  procedure  for  analyzing  multichannel  data  is  to  apply  Ozturk’s  algorithm  to 
each  channel  sq>arately.  The  capability  to  split  multichannel  data  into  its  component  channels 
already  existed  at  the  beginning  of  this  effort,  and  was  not  changed  here. 

In  iqjplying  Oztuik’s  goodness-of*fit  test  to  each  channel  of  a  multichannel  dataset,  the 
user  must  take  into  account  how  the  significance  level  is  altered.  For  example,  suppose  the 
data  ctmsists  of  two  stochastically  independent  channels.  Consider  testing  the  null  hypothesis 

■  Ho  :  xi(n)  is  distributed  Fi  and  X2(tt )  is  distributed  F2 
against  the  alternative 

•  Hi:  All  other  alternatives. 

Let  a  denote  the  significance  level  of  the  overall  test: 

a»Pr(Reject  Ho^Hotrue).  (3-111) 

This  statistical  test  can  be  considered  as  a  combination  of  two  independent  statistical 
tests.  The  hypotheses  fw  the  two  subtests  are 

•  Hi  :  xj(.n)is  distributed  Fj 

•  Hi:  Xj(n)  is  not  distributed  Fj , 

where  j  =  1,  2.  The  significance  level  for  these  subtests  is  a/ 

aj  ■  Pr  {Reject  Hi  I  Hi  true  ) .  (3-112) 

A  natural  decision  rule  for  the  overall  test  is  to  accq>t  the  null  hypothesis  if  both  subtests 
accq>t  tiieir  respective  nulls.  This  decision  rule  implies  Equation  3-113  for  relating  the  overall 
significance  level  to  that  of  the  two  subtests: 

1  ”  a  ^Pr  {{Accept  Ho  )  tad  {Accept  //^ )  I //q  true  )■  (1  -  ai)(l  •  02).  (3-113) 

Mote  generally,  if  a  statistical  test  is  the  result  of  J  indq>endent  subtests  with  identical 
significance  levels,  the  significance  level  of  the  subtests,  as  a  function  of  the  overall 
significance  level,  is  given  by  Equation  3-114: 

-I-(l-a)«'.  (3-114) 

This  relationshq)  is  not  encoded  in  the  MSPSS.  Rather,  the  user  must  account  for  it  when 
entering  desired  significance  levels.  The  user  must  also  account  for  any  correlation  between 
channels  and  the  conversion  of  complex  data  to  real  data. 
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3.4Jt  Real  Components  of  Complex  Data 

The  program  implementing  Ozturk’s  algorithm  accepts  single  channel  or  univariate  data 
created  by  one  of  the  other  programs  in  the  MSPSS.  This  data  can  be  teal  or  complex.  The 
user  is  provided  four  optimis  for  converting  this  data  if  it  is  cmnplex.  The  user  can  analyze 
the 

o  Real  part 

o  Imaginary  part 

•  Magnitude 

•  Angle. 

3.4  J  The  Quadratic  Form  q 

A  more  sophisticated  method  of  converting  multichannel  cmnplex  data  to  a  univariate 
quantity  was  al^  implemented  in  the  MSPSS.  This  method  is  based  on  the  theory  of 
Spheri^y  Invariant  Random  Processes  (SIRPs).  In  particular,  a  certain  quadratic  form  is  a 
sufficient  statistic  fn*  SIRPs,  and  Ozturk’s  algorithm  can  be  applied  to  many  realizations  of 
that  quadratic  form.  This  quadratic  fontn  is  further  defined  in  this  section. 

Let  X  ( 1),  X  (2).  ....  X  (IV )  be  a  multichannel,  possibly  cmnplex,  process.  Each  time  sample 
X  (a  )  is  itself  a  column  vector.  If  there  are  J  channels,  the  components  of  each  time  sample 
are  as  shown  in  Equation  3-115: 

xi(a) 
xt(n) 

(3-115) 

xj{n) 

The  entire  stochastic  [vocess  can  be  represented  as  a  single  vector  with  elements: 

x(l) 
x(2) 

(3-116) 

x(/V). 

The  JtJ  block  correladmi  matrix  summarizes  the  correlatitm  between  channels  and 
among  time  lags.  The  correlatitm  matrix  for  the  1th  lag  is  defined  as  in  Equation  3-117: 

^(/)-£[x(n)x"(n -/)].  (3-117) 

The  block  correlatitm  matrix  is  defined  by  Equation  3-118: 
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^(0)  Rd) 
R(.-l)  R(Q) 


R(N) 
R(N  -1) 


(3-118) 


R  m 


IRi-N)  R(l~N) 


R(0) 


The  quadratic  form  which  summarizes  a  SIRP  can  now  be  defined.  It  is  given  by 
Equation  3-119: 

qmx^R-^x.  (3-119) 

This  quadratic  form  is  always  real  and  positive.  Given  realizations  of  a  multichannel  radar 

return.  Oztuik’s  algorithm  can  be  applied  to  the  sequence  of  real  scalars  q\  q^ . q'**.  This 

data  has  the  form  needed  to  tq)ply  Ozturic’s  algoridim. 

The  quadratic  form  is  difficult  to  compute  in  practice.  Radar  returns  in  the  MSPSS  can 
mutain  iq)  to  4  chantielg  and  Up  to  20,000  time  sanq>les.  Thus,  the  matrix  R  can  be  up  to 
80,000  X  80,000.  (3nly  a  special  case  of  the  quadratic  form  was  implemented  so  as  to  reduce 
these  storage  requirements.  This  special  case  omsists  of  data  uncorrelated  across  time,  but 
possibly  correlated  across  channels. 

For  diis  special  case,  the  block  conrelatimi  matrix  R  becomes: 


(3-120) 


where  £  is  the  Jxf  covariance  matrix  showing  cross-channel  correlation.  The  quadratic  form 
is  dien: 


9 


£  Jt^(ii)£-‘jt(«). 


(3-121) 


3.4.4  Oztivk’s  Goodness  of  Test 

Let  X(i)  ^  s:  ...  ^  X(jv)  be  an  ordered  random  sanq>le.  These  values  can  be  sorted 
quadratic  forms  as  described  above,  a  single  chaimel  oi  rnd  data,  die  real  part  of  a  single 
rfiorniftl  of  complex  data,  etc.  The  applicatitm  of  Ozturk’s  algoridim  as  a  goodimss  of  fit  test 
allows  the  user  to  test  whether  or  not  diis  sample  is  firom  a  specified  distribution.  That  is, 
Qzturk’s  algorithm  provides  a  statistical  test  for  deciding  between  the  null  hypodiesis 

•  H^:  Xn  a  firom  a  qiedfied  fx'obability  distribudmi 
and  the  alternative 

•  Hi’,  otherwise. 

The  null  hypodiesis  need  (mly  be  specified  as  far  as  the  type  of  the  distribudtm  and  shape 
parameters.  Location  and  scale  parameters  can  be  any  values.  Normalization  removes  any 
effect  from  location  and  scale  parameters.  The  null  hypothesis  distribution  should  reflect 
whatever  conversions  were  used  to  gotierate  the  random  sample  to  which  Ozturk’s  algorithm 
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is  applied.  For  example,  suppose  the  data  being  tested  is  the  magnitude  of  a  complex 
Gaussian  process.  The  null  di^budon  shouM  be  Rayleigh  (WeibuU  with  a  shi^  parameter 
of  two),  the  magnihKle  of  a  complex  Gaussian  random  variable  is  Rayleigh. 

Ozturfc’s  statistic  Qn  consists  of  the  ordered  pair  (l/jy,  V^)  defined  by  Equatitms  3-122 
through  3-123: 


"  1  «i 

(3-122) 

"  •  -1 

(3-123) 

0,  -  A «!>(«(,)). 

(3-124) 

where  y(,)  are  the  standardized  values  for  dm  ordered  sample: 

Z(.)-X 

y<i)  •  , 

(3-125) 

X  and  i,  are  the  sample  mean  and  sample  standard  deviation  of  the  randmn  sample  and  is 
the  eiq)ectBd  value  of  die  ith  order  statistic  from  the  reference  distribution;  that  is,  the 
standard  Gaussian  distribution^.  is  the  Cumulative  Distribution  Function  for  the 

reference  distribution,  as  given  by  Equatitm  3-126: 

(3-126) 

Ozturk’s  algorithm  can  be  visualized  as  tracing  a  path  out  in  a  two  dimensional  space 
{U,  V),  as  shown  in  Hgure  3-7.  The  data  defines  die  magnitudes  ly(a)l  of  a  sequence  of 
-vectors  in  this  qiace,  vdiich  are  then  standardized  by  the  factor  {UN).  The  sequeime  of 
angles  between  these  vectors  ate  defined  by  the  Gaussian  reference  distribudon.  The  ordered 
pair  given  by  Ozturk’s  algorithm  is  die  endpoint  of  die  path  formed  by  linking  die  vectors 
together.  The  null  hypothesis  defines  the  expected  ent^ioint  of  this  distributitm,  as  well  as  the 
statistical  variation  a^t  diis  em^xnnt  and  ^  statistical  variation  of  the  padi  to  this  endpoint 

The  expected  order  statistics  from  the  Gaussian  reference  distribution  are  found  by  a 
Monte  Carlo  esdmatitm  algotidun.  Let  wq)  i  W(3)  ^  ...  ^  be  an  ordered  sample  from  a 
Gaussian  distribution  with  zero  mean  ud  unit  variance.  A  user-qiecifed  number  of 
realizationa  of  these  order  statistics  are  syndiesized.  Estimates  of  die  eiqiected  values  d  the 
enter  statistics  are  found  by  averaging  over  these  realizations  (see  Appendix  D). 

The  user  must  know  die  distribution  of  Qn  under  the  null  hypothesis  to  decide  between 
the  null  and  alternative  hypotheses  at  a  specified  significance  le^.  This  distribution  is  not 
known  in  closed  form.  Tte  expected  value  of  the  Ozturk  statistic  is  also  found  by  Monte 
Carlo  simulation.  Let  Z(|,  ^  ^  ...  ^  Z(W)  be  an  ordered  sample  fixim  the  null  hypothesis 

distrihutitm.  Many  realizations  of  tins  ordered  sample  are  generated.  The  magnitu^  of  the 
linked  vectors  for  each  realization  are  given  by  Equation  3-127: 


'a  omm  (iMnl  vanioa  of  Oawfc*i  alforidan  woold  allow  dia  umt  io  i^pecUy  iht  mfeNMa  diMribotian  m  w«U  at  tka  null  diaifibu- 
Soo.  Oonric  hna  abowa.  lionwver.  Hut  the  nomial  dtalribullan  mipiriaBy  worka  bati. 
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(3-127) 


The  sample  mean  and  variance  are  estimated  anew  for  each  realization.  The  angles  between 
the  vectors  do  not  change;  they  are  based  on  the  expected  values  of  the  ordered  sample 
fh»n  the  reference  disffibution,  as  described  above.  The  expected  value  of  Ozturk’s  statistic 
is  found  by  averaging  over  many  realizations  of  the  endpoint  of  these  linked  vectors.  The 
aven^  path  is  shown  in  Figure  3-7. 

Figure  3-7  also  shows  confideiue  regiois  for  Ozturk’s  statistic.  These  are  needed  for  a 
formal  statistical  test  The  confidence  re^tm  should  be  determiiMd  for  a  level  1  -  a,  where  a 
is  the  significance  level  of  the  test  If  Ozmrk’s  statistic  for  the  data  lies  outside  die  confidence 
region,  the  user  should  reject  the  null  hypothesis  that  the  hypothesized  distribution  describes 
the  data,  hi  the  example,  the  data  is  rejected  as  coming  from  a  Gaussian  distiibuticHi  at  the 
10%  level,  but  accqit^  at  the  1%  level. 


-  Data  -  Mean 

-  90%  -  99% 


Figure  3-7:  Ozturk’s  Algorithm 
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CalcuUdoa  of  confidence  regions  is  a  fairly  complex  procedure.  Two  methods  are 
stqjplied.  One  is  based  on  the  assumption  diat  the  ordered  pair  {Un,  Vm)  is  approximately 
bivariate  normally  distributed.  The  otbtf  is  baaed  on  transfonnations  from  the  Johnson  family 
of  distributions.  These  procedures  sre  explsined  in  mote  detail  in  Appaidix  F. 

3.4.5  Oiturk*s  Algortthm  as  ParametM'  Estimation 

The  Oztuik  algorithm  traces  a  path  in  a  two  dimensional  space  whose  endpmnt  is  the 
statistic.  Bach  distrilnitioa,  vdien  its  shape  parameters  are  specified,  if  any,  results  in  a  definite 
ejq)ected  enc^joint  These  ent^xxnts  can  te  compared  with  Ozturk’s  statistic  as  determined 
firom  an  obao^red  random  sample.  Under  this  theory,  the  distribution  that  best  fits  observed 
data  can  be  determined  baaed  on  Ozturk’s  statistic. 

The  statistic  is  calculated  from  the  data  and  cotqpared  with  the  endpoints  from  various 
possible  distributions.  The  endpoint  closest  to  the  observed  statistic  is  identified,  and  the 
distribution  associated  with  tiiis  endpoint  is  selected  as  the  best  fit  to  the  data.  Figures  3-8, 
3-9,  and  3-10  show  a  "distributimi  ^^>roximati(m  chart"  needed  for  this  purpose.  In  additimi 
to  the  patii  of  linked  vectors  determined  by  the  data,  the  esqjected  endp^ts  are  plotted  for 
several  distributions.  A  single  point  is  plo^  when  the  Prol^ility  De^ty  Functitm  for  the 
distributitm  omtaina  no  shape  parameters,  (»ly  a  location  or  scale  parameter  (Figure  3-8). 
Distributions  with  a  sii^le  shi^  parameter  result  in  a  curve,  approximated  by  a  finite  number 
of  points  along  die  curve  (Figure  3-9).  The  K  disdibution  in  this  figure  is  for  a  real  process, 
while  the  Gumbel  distribution  is  of  Type  n.  Distributions  with  two  shape  parameters,  such  as 
die  Beta  and  Johnstm  SU  distributions  in  Figure  3-10,  are  shown  by  a  frutdly  of  curves.  Each 
curve  fx  the  Beta  distribution,  fx  instance,  corresponds  to  a  different  value  of  one  shqie 
parameter.  The  odier  shape  parametx  varies  along  die  curve. 

The  distance  between  the  data  and  the  eiqiected  value  of  die  statistic  fx  a  given 
distribution  and  slupe  parameters  is  the  Euclidean  distance.  A  more  precise  method  would  be 
based  on  probabilities.  The  Euclidean  norm  is  reasonable  if  most  confidence  r^ons  are 
qiproximately  circular.  The  closest  distribution  is  then  diosen  as  die  desired  distribution  fx 
qiproximating  the  observed  data. 

3v45.1  Kstnaatiag  Shape  ftmcters 

Slupe  parameters  fx  die  selected  distribution  can  be  determined  frxn  Ozturk’s 
qproadmation  chart  A  distribution  with  xie  slupe  parametx  will  be  rqxesented  by  a  single 
curve  in  die  qproxunation  chart  The  shape  paianrelx  can  be  estimated  fix  such  a 
distribution  by  finding  the  point  along  this  curve  closest  to  the  value  of  Ozturfc’s  statistic  as 
calculated  from  the  data.  In  practice,  diis  curve  will  nX  be  known  exacdy.  Instead,  a  finite 
series  of  points  akmg  the  curve  will  be  known.  Thus,  the  shipe  parametx  estimate  must  be 
qproximiited. 

Sippose  (1^1,  Vi)  and  (I/j,  Vj)  in  Hgure  3-11  are  points  along  a  curve  corresponding  to  a 
distribution  with  one  shape  parametx  y.  Let  (U,  V)  be  a  value  of  Ozturk’s  statistic  calculated 
firom  the  data.  The  distribution  curve  is  ^^oximated  by  a  straight  line  between  {Uu  Vi)  and 
(Ub  V}).  The  problem  is  to  {X'oject  ((/,  V)  onto  the  line  at  point  (l/o,  Vo)  and  determine  the 
cxrespxiding  shape  parametx  fx  this  point. 
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OzturUc's  Stat'st‘c 


-0.2  -0.1  0.0  0.1  0,2 
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Figure  3-8:  Location  and  Scale  Distributions 
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Ozturk's  StatTstfc 
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Figure  3-9:  Distributions  with  One  Shape  Parameter 
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Ozturk'  s  Stat  T  st  ic 


Figure  3-10:  Distributions  with  Two  Shape  Parameters 
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(t/.V) 


Figure  3-11:  Linear  Interpolation  for  Shape  Parameter 
The  slope  of  the  line  between  (Ui,  V,)  and  (1/2,  V2)  is  m: 


£/a-£/, 

(3-128) 

Toe  equations  for  the  two  Unes  whose  intersectitm  is  (Uo,  Vo)  are: 

y  mmx  +(V,  -m  t/,) 

(3-129) 

y  x+(V  +  i-(/) 

(3-130) 

One  can  substitute  (Uo,  into  these  two  equations.  The  solution  of  the  resulting  system  (tf 

equations  is: 

(3-131) 

m^V  *m(U 

*  1  « 
m*+  1 

(3-132) 

An  estimate  of  the  conesponding  shape  parameter  is  then  found  by  iin».yr 
between  tiie  em^Kunts: 
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t/j -  i/o  f  Vi-Ut,' 


Yz 


(3-133) 


where  yi  and  are  the  shf^ze  parameters  conresponding  to  the  endpoints. 

Estimating  distributions  with  two  shape  parameters  is  more  complicated.  One  method 
would  be  to  use  linear  interpolation  based  on  three  points  around  the  obsoved  statistic.  The 
method  implemented  in  the  MSPSS  is  simpler.  Each  curve  is  treated  as  a  sqiarate  distribution. 
The  curve  closest  to  the  observed  statistic  is  found.  This  determines  one  sluq>e  parameter.  The 
other  is  found  by  linear  interpolation  along  the  selected  curve,  as  above. 


3A5,2  Eatiwating  Location  and  Scale  Fanuneters 

Ozturk’s  distribution  approximation  chart  can  be  used  to  identify  the  distribution  that 
most  closely  matches  an  oteerved  random  sample.  It  can  also  be  used  to  determine  shape 
parameters.  Estimates  of  location  and  scale  parameters,  however,  are  not  available  from  the 
q>proximation  chart.  The  approach  describe  in  (Shah  93)  was  implemented  to  estimate 
location  and  scale  parameters  in  a  manner  consistent  with  Ozturk’s  algcxithm. 

Let  X(i)  ^  X(2)  ^  ...  ^  be  the  ordered  random  sample  which  is  to  be  used  to  estimate 
location  and  scale  parameters.  Suppose  the  underlying  Probability  Density  Function  is 
/(x.ouP),  where  a  is  the  location  parameter  and  p  is  the  scale  parameter.  Consider  the 
standardized  order  statistics  i  =  1,  2, ...,  n: 

(3*134) 

P 

Let  p,  be  die  expected  values  of  these  standardized  (xder  statistics: 

M.-£IS(.)1  (3*135) 

In  practice,  p,  ,  i  =  1,  2,  ...,  n  is  found  by  Monte  Carlo  simulation.  The  order  statistics  S^y  are 
simulated  from  the  given  PDF  with  a  location  parameter  of  zero  and  a  scale  parameter  of 
unity  The  eiqzected  value  of  is  related  to  the  location  and  scale  parameters  and  the 
e}q)ected  values  of  S(,)  as  follows: 

£[X,„]-a  +  p.p.  (3-136) 

Consider  the  statistics  given  by  Equations  3-137  and  3-138: 

r,- £Coj(0i)X(,)  (3-137) 

i  >1 

T^•isin(^i)X^yy,  (3-138) 

I  - 1 

where  0,  is  as  defined  in  Equation  3-124.  These  statistics  closely  resemble  Ozturk’s  statistic. 
Their  expected  values  are: 

£[r,]-oo  +  ftp  (3-139) 


where 


£  l^jJ  -  c  a  +  4  P. 


(3-140) 
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(3-141) 


a  -  Cos  (0, ) 

I  >1 

6-tM.Co*(ei)  (3-142) 

■  -I 

cm  t,  Sin  (0. )  (3-143) 

i  -  1 

4- £|i,Sin(0i).  (3-144) 

imt 

It’s  fairly  ai^arent  by  symmetry  that  a  must  be  equal  to  zero.  Hence  estimates  of 
location  and  scale  parameters  are: 

p — r- 

,  EIT-jl-rfp 

a  • 

c 

The  expected  values  of  Ti  and  Tj  are  found  by  Equations  3-137  and  3-138  from  the  observed 
data.  As  has  already  been  explained  ii,  and  0;  are  found  by  Mtmte-Gtflo  simulation.  With  the 
determinatitm  of  the  location  and  scale  parameters,  the  Ozturk  algorithm  for  parameter 
estimation  is  completely  specified. 


(3-145) 

(3-146) 
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4.  IMPLEMENTATION 

Section  3  describes  new  capabilities  implemented  under  this  eif(»t.  These  capabilities 
were  implemented  with  several  interacting  programs  under  the  menu-based  subsystem  and 
new  analysis  sequences  under  the  UFI-based  system.  Each  program  under  the  menu-based 
system  corresponds  to  a  menu  choice.  Two  menu-based  programs  and  four  UFI-based 
analysis  sequences  were  designed  to  implement  the  senstM*  fusion  approach.  Kalman  filtering 
was  implemented  with  four  {U'ograms  and  four  analysis  sequences.  The  Extreme  Value 
Theory  resulted  in  the  modification  of  loglilrelihood  calculation  prt^ams  and  all  detection 
analysis  sequences.  The  Ozturk  algorithm  was  implemented  as  one  program. 

4.1  Sensor  Fusion 

The  senstv  fusion  t^roach  requires  the  following  capabilities; 

•  Synthesis  of  multicharmel  data  as  ARMA  processes  where  off-diagonal  elements  of  AR 
coefficient  matrices  are  zero 

•  Estimation  of  single  channel  AR  parameter  along  each  chaimel  of  multicharmel  data 

•  Estimation  of  the  error  covariance  matrix  after  tapped  delay  line  structured  prediction 
erTtH*  filters  are  q)plied  to  each  chaimel  of  multichamel  input 

•  A  sensor  fusion  filter  which  ^plies  a  tiyjped  delay  filter  to  each  chaimei  sq)arately  and  a 
second  phase  at  the  fusion  center  for  removing  interchaimel  correlation. 

Some  of  these  «q>abiiities  can  be  viewed  as  special  cases  of  functions  that  existed  at  the  start 
of  this  effort  In  particular,  the  previously  existing  tapped  delay  line  filter  has  parameters 
corresptmding  to  an  AR  model  of  its  input.  This  filter  is  fimctionaUy-equivalent  to  the  sensm* 
fusion  filter  when  the  AR  coefficients  are  diagonal: 

dll/)  0  0  ...  0 

0  did)  0  0 

0  0  did) 

A^d)m  .  (4-1) 

0 

0  0  0  ...  did) 

where  /  is  the  number  of  channels.  The  order  of  the  vector  AR  process,  where  each  single 
channel  is  an  AR  process,  is  the  maximum  of  the  mder  of  the  single  channel  AR  processes. 

To  take  advantage  of  existing  multichannel  processing,  all  files  of  AR  coefficients 
created  by  new  programs  were  stmed  as  diagonal  matrices.  This  design  decision  allowed 
previous  routines  for  synthesizing  ARMA  processes,  filtering  multichannel  inputs,  displaying 
AR  coefficients,  etc.  to  be  available  for  sensor  fusion  analysis. 

4.1.1  Estimation  Program 

A  menu-based  program  was  written  in  which  AR  model  parameters  are  estimated  using 
single  channel  estimators  for  each  channel  in  a  multichannel  process.  The  input  consists  of 
many  realizations  of  a  multichannel  signal  as  generated  by  the  multichannel  synthesis 
routines.  The  output  is  a  set  of  multichannel  AR  parameters  in  which  all  off-diagonal 
elements  are  zero.  These  parameters,  or  their  average,  can  be  stored  in  a  file  suitable  for 
reading  by  the  multichannel  tapped  delay  line  filter  and  other  programs  in  the  system.  In  other 
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words,  the  out{>ut  file  is  the  same  fonnat  as  the  output  file  for  the  multichannel  Yule-Walker 
estimatioa  program,  for  example. 

entry,  die  user  is  asked  for  the  name  of  the  ii^ut  file.  The  user  is  also  asked  to  input 
an  integer  vector  in  which  each  element  rqiresents  the  order  of  the  desired  AR  model  for  the 
correqxmding  chaniwl.  The  user  is  asked  to  choose  an  algorithm  for  estimating  the  AR 
parameters.  Cunent  opdmis  consist  of  the  Yule-Walker  and  Burg  algorithms,  these  being  the 
available  algorithms  already  implemented  for  single  channel  AR  estimation.  If  the  user  wants 
to  use  the  Yule-WaUcer  algmitfam,  the  user  must  indicate  whether  biased  or  unbiased  estimates 
should  be  found. 

For  each  realization,  the  program  estimates  AR  parameters,  including  the  variance,  fa*  an 
AR  model  for  each  channel.  These  estimates  are  stored  in  a  rqiresentatimi  of  a  multichannel 
AR  process  with  diagonal  matrices  for  AR  coefficients  and  the  covariance  matrix. 

The  method  of  displaying  these  estimates  for  each  realization  is  modeled  after  the 
already  existing  programs  for  estimating  multichannel  AR  coefficients.  The  user  is  given  the 
optitm  of  saving  the  estimates.  <»'  their  averages  over  many  realizations  to  a  specified  file. 
Variances  are  calculated  and  displayed  for  these  averages.  No  option  exists  for  calculating 
error  variances  based  on  the  actual  values.  Figures  4-1  and  4-2  show  an  example  of  user 
interaction  with  this  program. 

4.U  Covariance  Matrix  Estimation 

A  menu-based  program  was  also  written  to  estimate  the  covariance  matrix.  The  input  is  a 
multichannel  signal.  This  input  could  be  the  output  of  a  tapped  delay  line  structured  predictitm 
error  filter  in  which  the  AR  coefficients  are  as  determined  by  the  AR  estimation  program.  The 
matrix  used  to  remove  interdaaimel  correlation  in  the  second  phase  of  the  filter  should  be  set 
to  identity.  By  using  the  menu-based  system  in  this  fashiem,  tite  structure  of  the  sensor  fusion 
qiproach  is  preserved.  But  the  covariance  matrix  estimation  program  could  be  used  in  other 
instances  as  well. 

The  program  estimates  the  covariance  matrix  for  each  realization  of  the  input  process, 
e.g.  [xediction  error  filter  ou^ut.  These  estimates  and  their  averages  across  trials  are 
di^layed  to  the  user.  The  estimates  are  calculated  as  in  Equation  3-27  and  3-28.  The  user  is 
given  the  option  of  saving  the  average  estimates  to  a  file.  If  the  user  chooses  this  optitm,  the 
user  is  given  the  option  of  storing  these  covariance  estimates  with  AR  parameters  read  from 
an  already  existing  file.  This  furthtf  optimi  can  produce  a  coefficient  file  suitably  formatted  to 
serve  as  the  parameters  for  the  sensor  fusion  filter  shown  in  Figure  3-1.  Figure  4-3  shows  an 
exanqile  execution  of  the  covariance  estimatim  program. 

4.1.3  Sensor  Fusion  Analysis  Sequences 

Four  new  analysis  sequences  were  written  to  implement  sensor  fusion  in  the  UFI-based 
subsystem.  Figure  4-4  shows  the  UFI  menu  fw  sensor  fusicxi  sequences.  The  covariance 
matrix  estimation  sequence  merely  generates  a  user-specified  number  of  trials  of  white  noise 
and  estimates  the  covariance  matrix  from  those  realizations. 

The  sensor  fusirni  AR  parameter  estimation  sequence  is  more  complicated.  A  first  set  of 
realizations  of  an  AR  process  are  generated.  These  realizatitms  are  used  to  estimate  AR 
parameters,  where  each  channel  is  modeled  as  a  single  channel  AR  (xocess.  Once  the  average 
of  these  AR  parameters  is  determined,  a  new  set  of  realizations  of  the  AR  process  is 
generated.  Thew  realizations  are  filtered  by  a  tapped  delay  line.  Each  channel  is  filtered 
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MYWB  --  rhmmmi  AR  pinmeters  by  Ytile^Walker  or  Bntg  tlgoritfam. 

VenknU 

Iniiat  file  name  ?  riv 

Title:  AR.  a  « (-0.625.  -1.629).  (0l2S.  0.81) 

Output  fila  itHMt  7 

No  fib  name  entered,  is  dus  okay  ?  y 
Display  die  esthnates  on  the  tenninal  (y/h)  [Y]  ? 

Do  yon  want  to  use  the  Ynle-Walker  or  Bvg  method  (y/b)  7 
Enter  ”y"  to  ore  die  singie-channel  Ynb-Wafiger  algoritfam  to 
AR  preameteia  for  each  chaniid.  Eater  '’b'*  to  ore  the 
Bmgalgarid^ 

Do  yon  want  to  ore  the  Yule-Walker  or  Borg  method  (yfly)  7  b 
Enter  die  onbr  of  die  AR  process  for  each  channel  7  2, 2 
Trial  1  size:  20000 
Covariance  matrix  estimate: 

(1.06s9e«00.0.0000e400)  (O.OOOOb400.0.0000mOO) 

(0.0000e400.0.0000e400)  (3.2686e-Q2.0.0000MOO) 

AR  coefficient  estimates: 

A(l)-(-6J876e-01.4.9529e-03)  (0.0000e+00.0.0000e400) 
(O.OOOOe400.0.0000MOO)  (-1.6330e+00.-3J726e-03) 
Aa)-(14834e-01.-9.8383e-04)  (O.OOOOe+OO.O.OOOOMOO) 
(0.0000e+00.0.0000fr«00)  (8.1419e4)U.7286e-03) 

I¥oo^  Next  trial  All  trials,  or  Quit  (n/a/q)  [N1  7  a 
Trial  2  size:  20000 
Covariance  matrix  estimate: 

(l.Q574»»OO.O.OOOOe400)  (O.O000e4OO.O.00OOMO0) 

(O.OOQOmOO.O.OOOOmOO)  (32898e-Q2.0.0000e400) 

AR  coefficient  estimates: 

A(l)-(-6218Se4)1.4.2076e4)3)  (O.0000MOaO.0000»*00) 
(0.000Oe4OO.O.000OMO0)  (•Iri302et00.-13750e^)4) 
Aa>>a4799e4)lZ3717e4)3)  (0.0000e»00.0.0000e400) 
(0.000(b»O0j0.000OMO0)  (8.0997e4)U.30S6e4)4) 

Trial3sae:  20000 
Covariance  matrix  estimate: 

(1.0562ti400.0.0000e400)  (0.0000e400.0.0000e«00) 

(0.000(b«00,0.000(b»00)  (3  Js26e-az.0.0000e400) 

AR  t 

A(l)-(*6J2QSe4)i5.1982e4)3)  (O.OOOOmOO.O.OOOOmOO) 
(O.OOOOe*OO.O.OOOOefOQ)  (-1.6264et00.7J044e-04) 
Aa)-CLS694e4)l.-6.7751e-03)  (O.OOOOMOaO.OOOOe400) 

(aOOOOiMOO.O.OOOOe^  (8.0736»01.6.6720b4M) _ 

Figure  4-1:  Sensor  Fusion  AR  Estimation 
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i  nai4size:  zuuuu 
CovflriMDOB  outrix  csdouttB: 

(1.0624e4O0.O.0Q0OMOO)  (O.OOQOmOO.O.OOOOmOO) 

(aooooMoaaooooMO))  o  J659iB-az.o.ooooMOO) 

ARcoefficieatestiiiMies: 

A(lH-6J828e-01.-3.Q3Q5e-03)  (O.OOOOmOO.O.OOOOb400) 

(O.OOOOB400.0.0000BfOQ)  (•1.6307e+0ai.48S9)B-(») 

Aa)-a5S«9e-01.-3J2103e-04)  (O.OOOOMOaO.OOOOMOO) 

(aOOOOMOO.aOOOOD^  (8.1168e^l.-1.2603&O4) 

Trial  5  size:  20000 
Covariance  matrix  esttmaie: 

(1.0618D400.aOOOOe400)  (0.000(b400.0.0000MOO) 

(O.OOQOe400.0.0000B400)  (3^1e-a2.0.0000MOO) 

AR  ooefficieiit  esrimaies: 

A(lM-6J449e-01.-1367Se-Q2)  (O.OOOOM<X).O.OOOOe400) 

(0.0000e400.0.0000e»00)  M.6273ef00.83958e-04) 

Aa)-a4287e^l.7.6110e-03)  (O.OOOOMOO.O.OOOOa400) 

(O.OOOOe+OO.O.OOOOMOO)  (8.0730e^l.-2.3590e-03) 

Mean  values  for  S  trials 
Mean  covariance  matrix  estimate: 

(L0607e400.0.0QOOe400)  (0.0000e400.0.0000e^) 

(0.0000e400.0.0000e400)  (3:2622e-a2.0.0000e400) 

Variance  in  covariance  matrix  estimates: 

1^12e^  aOOOOetOO 
O.OQ0OetOO4.2464e4» 

Mean  AR  coefficient  estimates: 

A(lH-6^"09e-OU3064e-04)  (O.OOOOe+OO.O.OOOQe400) 

(0.00C  >OO.O.OOOOe+OQ)M.629SefOO.-lJ064e4)4) 

A(2)-aS0i6e4)l  J.8Q54e4M)  (O.OOOOMOaO.OOOQe400) 

(O.0000e4O0.O.000QMO0)  (8.10iaa4)1.4.8S27e4)4) 

Variance  in  AR  coeffiaent  estimates: 

A(l>i9.ti361e-05  O.OOOOmOO 
aOQOOe+00  1.1122e4)5 
Aa)-6.1818e-05  O.OOOOmOO 

aOOOO^tOO  1.3396e-05 _ 

Figure  4-2:  Sensor  Fusion  AR  Estimation  (Continued) 
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MVAR  ~  Eidiiiales  imiHiclmiiid  covsiaace  matrix. 

Veniao  13 
bpot  file  oamB  7  riv 

TiflB:  CmnplBx  Gausaum  noiae.  mean  >  0.  var  >  (2. 1/2) 

Save  cvcRefie  oovariaoce  rarimaiea  to  a  file  (yja)  fY]  7  ■ 

Deqtiay  die  eathnetea  on  dm  lenmnal  (yAO  [Y]  7 
Trial  1.  Sine:  20000  samplea. 

Cbamiell:  Memi-aS743e^.-l.a237e-Q2) 

C3miei2:  Meaii-(43771e-03337S8e4}3) 

Covariances:  a0144e^.O.OQOOMOO)  (-SJ891e4».2.2511e4») 
(-S3891e4)3.*2.2Slle-03)  (5.04fi8e4)1.0.0000e400) 
Frixesa  Next  trial  All  trials,  or  Quit  (n^/q)  [N]  7 
Trial  2,  Size:  20000  sanqries. 

Ommll:  Meaii-(-9.7988e-033.7998e-03) 

Cliannd2:  Meaim<6.9riS8e-03.7.7491e4)3) 

Covariances:  (1.9965e+00.0.0000e+00)  (43871e>03.-1.17S9e4)4) 
(43871e4)3.1.1759e-04)  (S.Q2SSe^l.O.000OB+OO) 
ftooesa  Next  trial  All  trials,  or  Quit  (n/a/q)  [N]  7 
Trial  3.  Size:  20000  sanities, 
awmel  1:  Meaii-(S.4823e4)3.-133S3e-02) 

Clininel2:  Meait-(2.4676e4)333S79e-03) 

Covariances:  (2.0082efOOX).OOOOMOO)  (-6.8590e4)3.7.9171e4)3) 
(•6.8590e4)3.-75171e-03)  (4.9942e-01.0.0000B400) 
Brocesa  Next  trial  All  trials,  or  Quit  (a/tjq)  [N]  7 
Trial  4.  Size:  20000  samples. 

Channel  1:  Mean<23110e^.l.6887e4») 
ainiimi2:  ltfeanK-3.49Sle4)3.13439p4)4) 

Covarimnes:  a01Q2efOOjO.OOOOe»00)  (•3J123e4)4.-33313e4») 
(•33123e4)433313e4B)  (43Sfi2e4)1.0.0000e400) 
lYocees  Next  trial  All  trials,  or  Quit  (nMq)  M  7 
Trial  5.  Size:  20000  samples. 

Chamiell:  Mean-ad013e4)3.-1.1193e4») 

Clinmel2:  Mea8K-2.663Se-03.-23687e-03) 

Covariances:  aQ029e«OO.O.OOOOetOO)  (S3705e4)3.73S92e4)3) 
(S37Q5e4)3.-73592e4»)  (S.0278e-01.0;0000e400) 
Aooess  Next  trial  All  trials,  or  Quit  (n/iiAi)  {N]  7 
Mean  valoes  for  5  trials 

Channel  1:  Mea»-(6.7404e4)4.-33440e^}  Variaace-8.4734e4)5 
Channel  2:  Mean-(13304e4)33.4077e4»)  Varianoe-33430e4)5 
Covariances:  aOOfiS»tOO.O.OOOOMa))  (•43434e4)4Z8557»03) 
(•43434e4)4.-2.8SS7e413)  (5.0105e^l.0.0000e400) 
Variances  in  covariance  esrimates: 

Variances:  4.7790e^  S.fi330e4)5 
_ S.6330b4)S  1.3004e-05 _ ‘ 


Figure  4-3:  Estimating  the  Covariance  Matrix 
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Select  a  sensor  fusion  sequence 

(S)© 

BBM  II  1  III  III  III 

) 

) 

BPiei  II  1  II 1  II  Hill  II  —1^^— 

i 

(^7^  f  Fusion  detection  analysis  in  clutter  and  noise 

) 

^  „■  . .  . 

Figure  4-4:  Sensor  Fusion  Analysis  Sequences 

sqMniely.  The  resulting  multichannel  ouq>ut  is  then  used  to  estimate  the  covariance  matrix 
for  the  signal  produced  after  intertemporal  conelatioa  is  removed.  The  final  outout  of  diis 
sequence  is  the  average  of  these  covariance  matrix  estimates. 

The  sequence  for  detecdon  analysis  in  noise  is  based  on  the  structure  shown  in  Hgure 
2>1.  The  Fq  filter  merely  decorrelates  its  iiqnits  across  channels;  it  does  not  contain  a  tapped 
delay  structure.  The  covariance  matrix  used  for  this  filter  is  estimated  based  on  a  numbtf  of 
realizations  of  white  noise.  The  F|  filter  has  the  structure  shown  in  Figure  3-1  and  its 
parameters  are  estimated  as  in  the  sensor  fusion  AR  parameter  estimation  algorithm.  The 
detection  analysis  (n’oceeds  by  first  determining  the  diresholds  for  given  false  alarm 
probabilities.  This  requires  the  synthesis  and  filtering  of  many  realizations  of  nmae.  The 
probability  (rf  detectiai  is  estimate  finxn  many  realizations  of  signal  and  noise.  This  whole 
process  ia  rqieated  a  user-specified  number  of  times  so  as  to  be  able  to  estimate  means  and 
variances  for  the  thresholds  and  probabilities  of  detection. 

The  sequence  for  detection  analysis  in  clutter  and  noise  is  close  in  structure.  In  this  case, 
both  filters  ^ve  the  structure  shown  in  Hgure  3-1,  with  corresponding  coefficient  estimation 
loops.  The  ouq>ut  of  this  sequence  is  also  estimates  of  thresholds  and  prdiabilities  of 
detectum  for  given  false  alarm  probabilities. 

4J^  Nfodel  Identification 

Four  programs  were  added  to  the  menu-based  subsystem  in  suf^xat  of  state  q>ace  model 
identification.  These  programs  convert  AR  model  parameters  to  state  space  model  parameters, 
synthesize  state  space  output  processes,  estimate  state  space  model  parameters  from  many 
realizations  of  the  output  process,  and  implement  a  Kalman  filter.  Four  sequences  were  added 
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to  the  UFI-baeed  subsystem. 

4J2.1  Coamskm  of  AR  Piaraineters  to  State  Space  Parameters 

The  program  for  converting  AR  parameters  to  state  space  parameters  allows  the  user  to 
either  specify  a  file  containing  Autor^resaive  model  parameters  or  enter  diem  direcdy  (Figure 
4-S).  These  parameters  are  the  lagged  wei^ts  in  the  AR  process  and  the  covariance  matrix 
for  the  driving  noise.  The  input  file,  if  any,  is  produced  by  one  of  the  other  programs  in  the 
system  such  as  the  AR  process  synthesis  or  Yul^Walker  estimation  programs. 

As  described  in  Section  3.2.I.3.  an  AR  iM‘ocess  can  be  modeled  as  an  innovations 

representati<m  of  a  state  ^ace  model.  The  program  calculates  die  innovadtms  model 

parameters  F,  K,  and  H"  as  described  there.  The  covariance  matrix  for  the  iimovations 

process  is  unchanged.  The  state  space  parameters  are  output  to  the  uso*  and  to  a  user- 
specified  file.  The  output  file  is  formatted  sudi  that  it  can  be  read  by  the  state  ^ace  jx-ocess 
syndiesis  or  Kalman  filter  programs. 


Gxivetti  AR  paiameiers  to  state  space  parameters. 

Venkm  1 J 

Read  AR  parameters  coefficients  from  a  file  (yM)  fYl  ? 
hqiot  file  name  ?  rivar 
Title:  Actnal  AR  coefficients 
AR  coefficients: 

A(l)-(-6J50Qe-01.0.0000e400)  (O.OOOOMOaO.OOOOe4€0) 

(O.OOOOmOO.O.OOOOmOO)  (-1.6290e40aO.OOOOMOO) 

A(2HZS000»O1.0.0000m00)  (O.OOOOMOO.O.OOOOMa)) 

(P.OOOOe400.0.0000MOQ)  (8.1000e4)1.0.0000e«00) 

Driving  *»«»»»  covariance  matrix: 

(l.QS4688afOO.O.OOOOOQe400)  (O.Q(XX)OOe400.0.00000Qe44)0) 

(aOOOOOOBfOO.O.OOOOOQe400)  (3.267Q52e-a2.aOOOOOQMOQ) 

Ate  dieae  panmeters  okay  (yM)  [Y]  ? 

Output  fito  name  for  state  qiace  pameter  ?  rlvaa 
Ude  of  thia  dataaet  7  State  Space  parameters 
System  Dynamics  matrix  F: 

(6JSaQOOB41.0.(IOOOOOe«00)  (O.OOOOOOe400.0.000(XnB400)  (-is00000e-01.0.000000»t()0)  (0.000 
(a000000ef00.0.000000e400)  (1.629000e^.O.000000B4OO)  (O.OOOOOOe«OO.O.OOOOOOMOO)  (-8.099 
(1.000000MO0.O.000000e4OO)  (O.000Q00e4O0A00000Q^»OO)  (0.00000Ob4O0.O.000000mO0)  (0.0000 
(aOOOOOO»tOO.O.OOOOOOMOO)  (1.00000OmOO.O.00000QmOO)  (0.(XX)OOQe400.0.000000»fOO)  (0.0000 

Hermitian  traaspoae  of  Observarion  matrix  H: 

(&2S0000»O1.0.Q00000e+O0)  (O.00000Oe4OO.O.00000QMOO)  (-2300000e-01.0.000000MOO)  (0.000 
(aOOQOOOMOO.O.OOOQOOe«00)  (1.62900OmO0.O.00000Q»4OO)  (0.000000&400.0.0QOOOOa*00)  (-8.099 

XalinMi  fiain* 

(1.000000MOO.O.000000e4O0)  (0.000(X»e400.0.000000B400) 

(a000000ef00.0.000000»»00)  (1.000000ei00.0.000000e400) 

(0.000000e400.0.000000e400)  (O.OOOOOOefOO.O.OOOOOOB^) 

(O.OOQOOOMOO.O.OOOOOOe400)  (0.000000e400.0.000000e400) 

Covariance  matrix  for  innovatioos: 

(1.0s4688e^.0.000000e400)  (0.000000e^.0.000000e400) 

(0.000000e400.0.000000w00)  (3.2670s2e-a2.0.000000w00) _ 

Figure  4-5:  Converting  AR  to  State  Space  Parameters 
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4 State  Space  Process  Synthesis 

The  state  space  process  synthesis  program  (Figure  4-6)  asks  the  user  to  either  enter  the 
model  paramettfs  dirMtly  or  provide  tlum  in  an  input  file.  These  parameters  omsist  of  F,  K, 
and  in  die  innovations  rq»esentation  of  the  state  space  model.  The  covariance  matrix  of 
the  innovadois  Q  must  be  specified  also.  The  state  q>ace  parameters  are  used  directly  in 
process  synthesis;  to  date  no  "shaping  fiinctitm"  ^iproach  has  been  defined  or  implemented. 


SyntliBsizes  stale  ^Moe  process 
Versian  1.4 

Do  you  want  to  set  die  pseudo  randooi  generator  seeds  (yM)  [N]  ? 

Enter  the  nntnber  of  tnt^  of  the  signal  to  be  generated  7  5 
Read  state  ^aoe  parameters  firam  a  file  (yM)  [Y]  7 
bqmt  parameter  name  7  rivas 
Title:  State  Space  parameters 

Deoompositiao  mtthod:  Cholesky.  LDU.  or  SVD  (c/l/s)  [C]  7 

Do  you  want  the  driving  noise  to  be  a  Gaussian  or  SIRP  process  (g/s)  [G]  7 

Lengdi  of  the  generated  signal  in  points  7  10000 

Enter  the  number  of  points  to  be  generated  faefom  saving  data  7  2000 

Driving  noise:  real  ody.  imaginary  only,  or  coioplex  (rfi/c)  7  r 

A  2  real  order  4  state  qiace  process  will  be  generated. 

System  Dynamics  matrix  F: 

(6js0000e-01.0.000000»t00)  (0.000000e400.0.000000e400)  (-2 JOOOOOe-01 .0.000000^00)  (0.000 
(aOOOOOOefOO.O.OOOQOOe400)  (1.629000e»00.0.000000e400)  (O.OOQOOOe+OO.O.OOOOOOefOO)  (-0.099 
(1.00Q000MOO.O.000000e4OO)  (O.OOOOOQe400.0.QOOOOOe400)  (O.OOOOOQe400.0.0QOQOO&»OQ)  (0.0000 
(a000000et00.0.000000e400)  (1.000000e400.0.000000e400)  (O.OOOOOQMOO.O.OpOOOOefOO)  (0.0000 

Hnrmitian  tranqxise  of  Observation  matrix  H: 

(6.2s0000e4)1.0.000000ef00)  (O.00000O&tOO.O.00000QMOO)  (-2j00000e-01.0.000000e+00)  (0.000 
(aOOOOOOMOO.O.OOOOOOe400)  (1.629QOO&tOO.O.OOOOOOe400)  (0.000000e400.0.000000e^  (-8.099 

iCaiman  Gain: 

(1.000000ef00.0.000000e400)  (O.000000e4O0.O.00000OMO0) 

(aOOOOOOe+OO.O.OQOOOOe400)  (LOOOOOOmOO.O.OOOOOOmOO) 

(a000000ie«00.0.000000b^)  (0.000000&»00.0.000000b^) 

(a000000ef00.0.000000e400)  (0.000000e400.0.000000e400) 

Covariance  matrix  for  itmovatians: 

(1.054688et00.0.000000e«00)  (O.OOOOOOB40aO.OOOOOOb400) 

(aOOOOOOe*OO.O.OOOOQO^»00)  (3267Q52e-a2.a000000e4OO) 

Choledcy  decompositioa  of  oovariance  matin: 

(l.Q26980MOO.O.OOOOOOe400)  (0.000000&»00.0.000000b400) 

(O.OOOOOOMOO.O.OOOOOOe400)  (1.8074996-01.0.000000^00) 

The  driving  noise  will  be  Gaussian. 

10000  time  sanqiles  will  be  generated,  with.  2(XX)  time  samples  priming  the  process. 

Are  these  parameters  okay  (y/n)  [Y]  7 
Onqxit  file  name  7  riv 

Title  of  this  dataset  7  Synthesixed  state  space  process  _ 

Figure  4-6:  State  Space  Synthesis 
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The  user  can  specify  the  covariance  matrix  directly  or  in  factored  form.  Factorization 
options  ccMisist  of  the  Cholesky,  L-D-U,  and  SVD.  The  program  does  not  provide  any  checks 
on  the  factors  the  user  enters.  For  example,  the  program  does  not  verify  whether  the 
eigenvector  matrix  provided  as  an  SVD  factor  is  Hdrmitian.  If  the  user  specifies  the 
covariance  matrix  directly,  the  uso*  is  asked  which  decon^rosition  should  be  used  in  process 
synthesis. 

The  user  can  choose  to  specify  Gaussian  or  K-distributed  Spherically  Invariant  Randtmi 
Process  (SIRP)  driving  noise.  The  user  must  also  q>ecify  the  random  number  seed,  the 
number  of  time  samples  in  a  single  realization  of  the  ou^ut  process,  the  numbo'  of 
realizations  (trisls)  of  the  process,  and  the  number  of  time  samples  to  first  generate  and  then 
discard  in  "priming"  each  realization.  This  last  input  allows  the  user  to  produce  steady  state 
output  larocesses  by  discarding  transient  behavior. 

The  program  generates  tfie  specified  realizations  of  the  process  and  saves  them  in  a 
user-specified  output  file.  The  output  file  can  be  read  by  many  other  programs  in  the  system, 
including  plotting,  statistics,  and  estimation  routines.  The  innovations  are  generated  as 

e(a)mT2  (n).  (4-2) 

where  z  (n )  is  either  Gaussian  or  SIRP  white  noise  uncmrelated  across  channels  with 
(diagonal)  covariance  matrix  D.  If  the  user  chooses  a  Cholesky  decomposition,  D  is  the 
identity  matrix  and  7  is  the  lower  diagonal  matrix  C  where 

dmCC"  .  (4-3) 

If  the  user  chooses  a  L-D-U  decmnposition.  f  is  the  lower  diagonal  matrix  L  with  ones  along 
the  main  diagonal  and  where 

dmLDL".  (4-4) 

If  the  user  selects  the  SVD,  D  is  the  matrix  of  eigenvalues  A  and  T  v&Q  where 

QmQAQ".  (4-5) 

If  the  user  entoa  the  model  parameters  directly,  the  user  is  given  the  option  of  specifying 
a  file  in  which  diey  are  saved.  This  option  permits  the  user  to  resynthesize  a  particular  case;  it 
is  also  useful  for  the  innovations  filter  program. 

4  J.3  Estimation  of  State  Space  Parameters 

In  die  program  (Figure  4-7)  for  estimating  the  parameters  of  the  state  space  model,  the 
user  is  asked  to  specify; 

•  The  iqiper  bound,  L.  on  the  size  of  the  block  state  space  vectm* 

•  Whether  biased  or  unbiased  estimates  of  conrelatimis  should  be  used 

•  The  model  order,  if  the  user  wants  to  input  it  diiecdy 

•  The  method  to  be  used  for  calculating  the  model  order 

•  The  bound  for  calculating  the  model  (xder,  if  the  user  wants  the  model  order  to  be 

internally  calculated 

•  Whedier  diagnostic  inframation  should  be  displayed  on  the  user’s  terminal 

•  The  file  containing  the  process  from  which  the  parameters  are  to  be  estimated 
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•  The  output  file  in  which  the  parameter  estimates  are  to  be  stored. 


SSC  -  state  space  model  pvameiers  nsiiig  the  Scientific  Studies  Cotporelioo  algaridnn 

VeniaalJ 

lopot  file  naote  7  rhr 

Title:  Syiaheauad  state  space  process 

Output  file  name  far  pawneter  estimates  ? 

No  file  name  entemd,  is  this  okay  ?  y 

Opper  bound.  L.  on  si»  of  state  space  Mock  vector  7  4 

Cahailate  biased  or  imbiaaed  lags  (u/b)  7  n 

MODEL  ORDBl  MENU 

0  •  Enter  the  model  onier  dnecdy 

1  •  Use  d»  canonical  cornel arions  to  calculate  d»  model  order 

2  •  Use  die  nnnnaliTnd  nmning  sum  of  the  canonical  coaeladons 

3  •  Use  the  squares  of  the  canonical  coneladaos 

4  -  Use  the  nonnalized  nmning  sum  of  squares 

5  •  Use  the  log  parameters 

6  •  Uae  the  mutual  infonnadoa  parameteR 

Model  order  cakalatian  method  7  0 

Model  order72 

Display  diagnostic  information  (y/n)  7  n 
Order  2 

System  Dynamics  matrix  F; 

a677S40e4)3.0.000QQOe»00)  (-3.43fi075e4)1.0.000000MOO) 

(5.0S4108e412.0.000000e^)  M.02919ie-01.0.000000B+OQ) 

Heraiitian  tranqpoae  of  Obsmvation  matrix  H: 

(&399S69iB4)2jaOOOOOOe»(n)  a:i90421e-Oi.aOQOQ(»B^ 

(1.4763S3e^O.OOOOOOMOO)  (-236663  le4)1.0.0QOOOOM«Q) 

teahii—  Gmn: 

(•3J26846e4>4.O.(XX)000e4(lO)  (13597020*00.0.0000000400) 

(1.319421o4)2.0.000000o400)  (13334S6e-01.aOOOOOOe40Q) 

Covatianoe  matrix  for  ianovatiflns: 

(1.16044So400.0.000000o400)  (1323281e-01.0.000000o400) 

(1323281e4)1.0.000000o400)  (4312954e-Ol.a00Q000»4O0) 


Figure  4-7:  State  Space  Model  Identification 


The  program  averages  the  estimates  of  die  lag^  correlatiaos  over  all  trials  to  generate  the 
needed  estimates  A,.  Using  this  infonnatitm,  the  program  estimates  the  state  space  model 
parameters  for  an  innovations  rqnesentation  using  the  algorithm  in  Sectitm  3.2.3. 
Specifically,  the  ivogram  generates  estimates  of  the  system  dynamics  matrix  F,  die  Hermitian 
tranqxMe  of  the  observation  matrix  H,  the  Kalman  gain  K,  and  the  innovations  covariance 
matrix  D.  The  model  order  and  these  estimates  are  echoed  to  the  user’s  terminal.  They  are 
also  saved  to  a  file  in  a  format  readable  by  the  iimovations  filter  (xogram. 

the  user  indicates  diagnostic  information  should  be  di^layed,  the  following  additional 
information  is  displayed: 
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•  The  singular  values  of  the  past  correlation  matrix,  Rp,l.l 

•  The  singular  values  of  the  future  correlaticm  matrix,  Rf.l.l 

•  The  correlations,  their  normalized  running  sum,  the  squares  of  the  canonical 

caneladoas,  tlw  normalized  running  sum  of  the  squares,  the  log  parameters,  and  the 
mutual  information  parameters  (Only  the  series  used  in  calculating  the  order  is  printed  if 
the  user  indicates  diagnostic  informatitm  should  not  be  printed) 

•  The  singular  values  of  the  system  dynamics  matrix  F 
4,2.4  InnovatHNis  Filter  Program 

In  die  Kalman  filter  implementadmi  (Figure  4-8),  the  user  is  required  to  qiecify  the  name 
of  a  file  containing  the  input.  This  file  contains  a  number  of  trials  of  the  process  in  the 
customary  format  for  multichannel  processes  in  this  system.  The  program  also  asks  the  user 
for  two  output  files.  One  file,  which  is  optional,  is  for  storing  the  estimates  of  the  state 
variables.  Hie  other  file  contains  the  decorrelated  innovations  filter  output. 

The  user  is  given  the  option  of  reading  the  filter  parameters  F,  K,  H",  and  Q  from  a  file 
or  inputting  them  directly.  Hie  user  is  asked  for  two  additional  parameters.  One  is  the  length 
of  dm  innovations  sequence  which  must  be  less  than  or  equal  to  number  of  time 
samples  in  the  input  process  N.  The  program  discards  the  first  N  -Ni  time  samples  of  the 
innovations  to  remove  transient  behavior. 

In  order  to  whiten  the  innovations  vector  across  channels,  the  covariance  matrix  a  is 
factored.  Factorization  options  omsist  of  Cholesky,  LDU,  and  Singular  Value 
Decompositions.  The  values  of  die.  resulting  channel  variances,  o/,  are  stored  in  a  user- 
supplied  file  for  later  use  in  calculating  a  Iqglikelihood  ratio  statistic. '  Given  the  user  ii^iuts, 
die  program  calculates  the  innovatimis  {x-ocess  by  die  Kalman  filter  as  described  in  Figure  3- 
6.  A  user-specified  number  of  time  sanqiles  of  the  innovatims  process  are  discarded.  The 
final  ouqwt  of  die  filter  is  decorrelated  across  channels  as  above. 

4,2J  Kalman  Filter  Analysis  Sequences 

Four  new  analysis  sequences  were  written  to  implement  the  state  space  model  and 
IGdman  filter  umler  the  UFI-based  subsystem.  Hgure  4-9  shows  the  UFI  menu  for  state  space 
sequences.  The  sequence  for  estimating  state  space  model  parameters  for  a  signal  generates  a 
number  of  realizations  of  an  AR  process.  Thew  realizatimis  are  used  to  estimate  state  ^ace 
model  parameters.  These  estimates  are  the  final  output  of  this  analysis  seqiwnce. 

The  sequence  for  estimating  state  space  modd  parameters  for  the  sum  of  signal  and 
noise  is  very  similar.  In  this  case,  however,  the  estimates  are  based  on  the  sum  ^  an  AR 
I^ocess  and  white  nmse. 

The  two  detecticm  analysis  sequences  closely  resemble  tme  anotho*.  They  implement  the 
signal  detection  algorithm  illustrate  in  Hgure  2-1.  Both  filters  Fq  and  Fi  are  Kalman  filters. 
The  parameters  of  each  (rf  these  filters  are  estimated  before  determining  thresholds,  given 
false  alarm  jx-obabilities,  or  probabilities  of  detection,  given  threshold,  ijke  all  other 
detection  analysis  sequences  in  the  UFI-based  subsystem,  this  sequence  encloses  the  whole 
process  in  an  outer  loop.  Thus,  the  final  ouq>ut  includes  means  and  variances  for  thresholds 
and  probabilities  of  detection. 
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MKAF  -  ^**«^»«**  ±e  imovatioiis  preoess  using  the  Kiimen  filler 

Venkn  1.4 

taiMt  file  name  7  riv 

71^'  Symhesoed  stile  spmx  process 

Rle  fat  deoonelaied  imovitiaas  oatput  ?  rhrl 

Titb  of  this  dataset  ?  Ineovatiaes 

Hk  for  oatpot  of  the  state  variables,  if  desired  ? 

No  file  name  entered,  is  this  okay  ?  y 

Hk  for  ooqiat  of  variances  of  the  deconelaied  innovations  7 

No  fik  nanre  entered,  is  this  okay  7  y 

Enter  the  index  of  the  first  ontpot  value  7  9971 

Do  yon  want  a  Qiolesky.  LDU.  or  SVD  (cfiM  covariance  decoopositian  [C]  7 
bpnt  the  Kalman  filler  parameters  directty  or  firam  a  fik  (d/D  [D]  7  f 
Parameier  fik  name  (Kalman  filire  mattkes)  7  rIvreZ 
Title:  Estimated  state  ^ace  parameteis 
A  2  Older  2  state  qwoe  process  will  be  filtered. 

System  Dynamics  matrix  F: 

a.d7754(k4B.0.000Q0(k^)  (•3.43607Se4)1.0.(XX)Q0(k+(X)) 
(S.Q54iO8fr4)2.O.(X)000Oe+O0)  (•l.Q29191e^l.0.(X)0(X)0M(») 

Hemitian  tnmqtose  of  (Dbservation  matrix  H: 

(8.399S69e-Q2.0.00000(kt(»)  (729042 le-01.0.(X)OOOOe+OQ) 
(1.4763S3e^O.OOOQOOefOO)  (-2^6663 le^l.O.OOFIQOOe^ 

gahiian  (jtin: 

(•3226846&<M.0.(XX)00Ob4O0)  (U597Q2e«OO.O.OQOOO(k400) 
(lJ19421e^O.OOQOO(kfOO)  (1.8334S6e^l.a0Q000(k+0Q> 

CovarianoB  matrix  for  innovations: 

(l.lfi(M4Sot00.0.0(IQQQ0e4<»)  (1  J23281e^l.aOOOOOOMO(9 
(lJ23281e4)1.0.000000»ta))  (4.S129S4e4)l.a000000Bt€Q) 

(Zboksky  decomposition  of  covariance  matrix: 

(1.077240frf00.0.00000(k400)  (0.00000(k400.0.00000(k400) 
(1228400e^l.0.000000e400)  (6.604S87e-01.0.000000e+00) 

Inverse  of  Choleiky 

(928298Se4)1.0.00QOO(k«QO)  (O.OOQOOOe4a}AOOO()0(k4()0!l 
(•1.7265fiOe^lj0.0Q00QQb4O0)  (1314Q99ie400.0.000(XXk400) 

The  first  9970  tiaie  ismplft  will  be  discarded. 

Length  of  innovations  setreence  was  30 _ 


Figure  4^:  Kalman  Filter 
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Figure  4-9:  State  Space  Analysis  Sequences 


43  Extreme  Value  Theory 

Implementing  die  Extreme  Value  Thetxy  (EVT)  did  not  require  writing  any  new  menu- 
based  programs  or  UFI  sequences.  Instead  the  programs  iwiplemftnting  loglikelikood  ratios 
were  modified,  along  with  the  sequences  that  perform  detection  analysis  with  "a  prioi" 
estimates.  Only  two  programs  calculate  loglikelihood  ratios  in  the  menu-based  system.  One 
implements  the  SIRP  and  Gaussian  multichannel  loglikelihood  ratio  statistics.  The  other 
implements  the  loglikelihood  ratio  for  the  unconstrained  quadrature  case.  The  detection 
analysis  sequences  were  modified  to  include  an  option  for  the  EVT. 

Three  preliminary  modifications  were  made  to  the  two  loglikelihood  ratio  programs  in 
the  menu-based  system  before  adding  the  EVT.  When  this  task  began,  these  programs  did  not 
calculate  a  threshold  fm*  a  given  false  alarm  probability.  Nor  did  they  a  detection 

fx'obability.  Rather,  in  the  threshold  estimation  mode,  the  user  was  asked  for  an  index  into  a 
sorted  list  of  thresholds.  In  the  detection  probability  estimation  mode,  the  (xograms  rqxxled 
the  number  of  loglikelihood  values  above  a  given  threshold  and  the  total  number  of 
loglikelihood  values  estimated.  Finally,  the  user  was  (xiginally  required  to  rerun  these 
programs  if  detection  probabilities  were  desired  for  different  thresholds,  even  if  these 
estimates  were  based  on  the  same  sample  of  loglikelihood  ratios. 

These  programs  now  report  a  detection  probability.  Let  X,  i  Xj  s  ...  s  X,  be  a  sample  of 
order  statistics  for  the  relevant  loglikelihood  ratio.  For  estimating  the  detection  probability, 
tlMse  estimates  should  be  generated  under  the  alternative  hypothesis  in  which  a  ^igrini  is 
present.  Let  q  be  the  user  defined  threshold.  Let  j  be  the  smallest  index  such  that  for  all 
k  ij,  Xj  >11.  Consequently,  n  -y  +  l  is  the  number  of  loglikelihood  values  strictly  greater 
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than  given  direshold.  Previously,  the  programs  rqnxted  n  -  j  i  and  r  to  the  user.  Now, 
diey  also  report  (n  -  j  +  Din,  which  is  an  estimate  of  the  probability  of  detection. 

The  programs  also  now  estimate  a  threshold  for  a  given  false  alarm  probability.  The 
inttipolation  procedure  incorporated  into  the  menu-based  programs  is  the  same  as  that 
previously  implemented  in  the  UFI-based  system.  In  this  case,  let  X|  ^  X2  ^  ...  ^  X,  be  a 
sample  ^  order  statistics  for  the  relevant  loglihelihood  ratio  generated  under  the  null 
hypothesis  without  a  signal.  Let  F  be  the  underlying  cumulative  probability  distribution,  and 
1^  p  denote  the  desired  false  alarm  {x-obability.  The  program  must  find  an  estimate,  such 
diat 

(4-6) 

where  is  an  estimate  of  11. 

Since  F  is  the  distribution  function  for  a  ctmtinuous  probability  distributitm, 

l-F(il)-Fr(X  >ii)-Fr(X  iti).  (4-7) 

Ignoring  ties,  the  number  of  sample  values  greater  than  or  equal  to  Xj^.i  is  n  -  j .  Thus,  for 
any  7  1,  a  reasonable  estimate  of  F  (X; .,.  1)  is 

F(X;^,)-1--^-^.  (4-8) 

n  n 

or, 

j^nFiXj^O-  (4-9) 

IfR(l-p)isan  integer,  Xj  ^  t  seems  a  reasonable  value  for  fi.  In  general,  this  value  cannot 
be  assured  of  being  an  integiex.  Therefore,  let  n  (l^-/> ) >7  ;  where  7  is  the  integer  part  of 
n(l-p)  and  0  ^  <  1.  The  estimate  of  the  threshold  is 

fi-(l-g)X>^,+«X^^2.  (4-10) 

This  is  the  estimate  currently  used  in  the  h^PSS.  The  above  argument  is  only  a  plausibility 
argument,  not  a  formal  derivation.  Equation  10.4  in  (Rangaswamy  93)  gives  a  sli^tly 
different  estimate  which  is  asymptotically  equivalent. 

The  two  loglikelihood  ratio  programs  now  have  an  additional  loop.  Previously,  if  the 
user  desired  to  calculate  thresholds  for  different  indices,  or  detectitm  probabilities  for  different 
thresholds,  the  user  was  required  to  rerun  the  program.  Now,  after  reporting  the  desired 
result,  tile  {X’ogram  asks  if  the  user  would  like  to  estimate  anotiier  threshold  or  estimate 
anotiier  det^on  probability,  whichever  is  ai^opriate.  This  optitm  is  rqieated  until  the  user 
indicates  no  more  estimates  are  desired.  Addititmal  estimates  are  found  based  on  the 
previously  calculated  sample  of  loglikelihood  ratios.  Since  most  of  the  computatitmal  time  is 
mqiended  in  these  programs  in  calculating  this  sample,  addititmal  estimates  are  reported  very 
quickly. 
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Appendix  A. 

Notation 

An  intennediate  matrix  in  the  Canonical  Correlations  Matrix. 

AR  weights;  each  AR  coefficient  is  a  matrix. 

Matrix  AR  coefficients  for  die  clutter. 

Ihe  Cholesky  decomposititm  of  £. 

The  Cholesky  decomposititm 

(1)  The  diagonal  matrix  in  the  LDU  decomposition  of  £.  (2)  A  matrix  in  the 
state  space  model. 

The  diagonal  matrix  in  the  LDU  decomposidtm  of  £« . 

The  system  dynamics  matrix  in  the  state  space  model. 

(1)  The  filter  corresponding  to  die  null  hypothesis.  (2)  A  Cumulative 
Distribudon  Funcdon. 

(1)  The  filter  corresptmding  to  the  altemadve  hypothesis.  (2)  A  Cumulative 
Distribudon  Funcdon. 

A  matrix  in  the  state  space  model. 

The  observation  matrix  in  the  state  q>ace  model. 

An  estimate  of  the  stochastic  block  Hankel  matrix  in  the  Canonical 
Correlaticms  Algorithm. 

The  column  shifted  block  Hankel  matrix  in  the  Canonical  Cmrelations 
Algorithm. 

The  null  hypothesis. 

The  alternative  hypodMsis. 

The  number  of  channels. 

The  Kalman  gain  matrix  in  the  innovations  rqiresentation  of  the  state  space 
model. 

The  anqilitude  matrix  for  the  clutter. 

(1)  The  lower  triangular  matrix  in  the  LDU  decomposition  of  L.  (2)  The 
number  of  lags  used  in  estimating  correlation  matrices  in  die  Canonical 
Correlations  Algorithm. 

The  lower  triangular  matrix  in  the  LDU  decompositimi  of  . 

The  dimension  of  the  state  vector  in  die  state  s^  odel. 

A  sample  mean  for  a  sample  size  of  k  items. 

The  number  of  time  samples  in  a  realization  of  a  stochastic  process. 

(1)  The  number  of  realizations  used  in  estimating  state  ^>ace  parameters  in 
die  Canonical  Correlations  Algorithm.  (2)  The  number  of  realizatimis  of  the 
loglikelihood  statistic  in  the  Extreme  Value  Theory. 
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Tbe  order  of  an  AR  iH’ocess. 

The  matrix  of  right  hand  eigenvectors  in  the  Singular  Value  Decomposititm 

ofL. 

Oztuiic’s  statiatic. 

The  matrix  of  ti^t  hand  eigenvectors  in  the  Singular  Value  Decompositioa 

ofE„. 

The  block  correlation  matrix  for  a  radar  return. 

The  correlation  matrix  showing  cross  channel  conelaticms  for  a  radar  return  at 
lag  /. 

The  correlation  matrix  showing  cross  channel  correlations  for  the  clutter  at 
lag  /. 

An  estimate  of  the  future  block  correlation  matrix  in  the  Canonical 
Correlations  Algorithm. 

An  estimate  of  the  past  block  correlation  matrix  in  the  Canonical  Correlations 
Algorithm. 

A  random  variable  used  in  generating  SIRPs. 

The  sample  variance. 

The  diagtmal  matrix  of  eigenvalues  in  the  Singular  Value  Decomposition  of 
the  matrix  i4  in  the  Cantmical  Correlations  Algorithm. 

A  square  matrix  in  the  Can<»ical  Correlations  Algorithm  formed  by  taking 
the  upper  M  rows  and  M  columns  of  5a. 

Ihe  diagonal  matrix  of  eigenvalues  in  tbe  Singular  Value  Decomposition  of 
the  estimate  of  the  future  blodc  correlatitm  matrix  in  die  Canonical 
Conelatitms  Algorithm. 

The  diagonal  matrix  of  eigenvalues  in  the  Singular  Value  Decomposition  of 
the  estimate  of  the  past  block  correlation  matrix  in  tbe  Canonical  Correlations 
Algorithm. 

The  sample  variance  for  a  sample  size  of  k  items. 

(1)  Either  C.  L,  or  Q  in  the  Cholesky.  LDU,  or  Singular  Value 
Decmnposition  of  £.  (2)  A  matrix  rqxesenting  a  basis  transformation  of  the 
state  space  in  the  innovations  re{»esentation  of  tbe  state  space  model. 

A  statistic  used  in  estimating  location  and  scale  parameters  in  tbe  Ozturk 
algorithm. 

A  statistic  used  in  estimating  location  and  scale  parameters  in  the  Ozturk 
algorithm.  ^ 

A  transformation  matrix  in  the  Canonical  CcMTelations  Algtxithm. 

A  transformation  matrix  in  the  Canonical  Correlations  Alg<xithm. 

The  sampling  period  for  the  clutter. 

Bither  C„,L^,or  Q„  in  the  corresponding  decomposition  of  . 
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The  matrix  of  right  hand  e^^vectors  in  the  Singular  Value  Decomposition 
of  the  matrix  >4  in  the  Qmooical  Correlatioas  Algorithm. 

A  component  of  Ozturk’s  statistic,  the  ordered  pair  (  Vh,  Vh  ). 

The  matrix  of  li^t  hand  eigenvectors  in  the  Singular  Value  Decomposition 
of  the  estimate  of  the  future  block  correlation  matrix  in  the  Cimonical 
Correlatioas  Algorithm. 

The  matrix  of  right  hand  eigenvectors  in  the  Singular  Value  Decomposition 
(tf  the  estimate  of  the  past  block  coneladm  matrix  in  the  Ouumical 
Correladoos  Algorithm. 

The  matrix  of  left  hand  eigenvectors  in  the  Singular  Value  Decomposidmi  of 
the  matrix  i4  in  the  Canonical  Correlatioas  Algorithm. 

A  conqxment  of  Ozturk's  statistic,  the  ordered  pair  (  Uh,  Vh  ). 

The  complex  conjugate  of  X . 

The  inverse  of  the  matrix  X. 

The  Hermidan  transpose  of  tlte  matrix  X. 

The  upper  square  matrix  in  the  Canonical  Correlations  Algorithm  consisting 
of  the  ^iprec^ly  ntm-zero  entries  in  the  matrix  ^ 

A  <Mie-element  lag  in  a  system  block  diagram. 

A  matrix  in  the  Canonical  Correlations  Algorithm  formed  from  the  first  J 
rows  and  M  columns  of  tire  matrix 

A  tail  value  of  the  ordered  loglikelihood  statistics. 

A  matrix  in  the  Canonical  Correlatioas  Algoridun  formed  from  die  first  M 
rows  and  J  columns  of  dw  matrix  7> Hix- 

A  scalar  AR  weight. 

An  estimate  of  an  AR  coe^dmit. 

A  parameter  used  in  generating  K  distributed  SIRPs. 

Multichannel  clutter  (e.g.  a  vector  AR  process). 

An  index  over  the  AR  wei^ts. 

The  lag  valtte  at  which  die  shaping  fimction  for  the  clutter  peaks. 

The  number  of  loglikelihood  statistics  in  the  tail. 

In  Ozturk’s  algorithm,  the  expected  valtte  of  an  order  statistic  from  the 
standard  reference  distributimi. 

(1)  Noise  (e.g.  Gaussian  white  nmse  correlated  across  channejs,  but  not  in 
time).  (2)  An  index  for  time  in  a  vector  stochastic  process.  (3)  The  sample 
size. 

A  false  alarm  probability. 

A  quadratic  form  that  is  a  sufficient  statistic  for  a  Spherically  Invariant 
Randmn  Process  (SIRP). 


70 


9' 

St  sin) 

«(«) 

win) 

X,  xin),  xin) 
x‘in) 

Sin  In  -  1) 

Jf. 

J 

xjin) 

yin) 

y<i) 

*(») 

A 

A 

A<‘> 

A, 


A, 

A, 

n 

£ 

£ 

i 


The  quadratic  form  q  as  calculated  for  the  ith  realization  of  the  stochastic 
process  modeling  a  SIRP. 

A  multichannel  signal  (e.g.  a  vector  AR  [vocess). 

A  sample  standard  deviadoo. 

The  iiq)ut  process  in  the  state  space  model. 

(1)  Noise  (e.g.  Gaussian  white  noise  correl^ed  across  channels,  but  not  in 
time).  (2)  ^asurement  imise  in  the  state  space  model. 

In  Ozturk’s  algorithm,  an  order  statistic  from  the  standard  reference 
distribution. 

A  multichannel  radar  return  (vector  stochastic  process). 

The  nth  time  sample  from  tlm  ith  realization  of  the  stochastic  process  xin). 
An  estimate  of  a  multichannel  radar  return. 

An  element  of  a  random  sample. 

An  element  of  a  ordered  randtnn  sample,  an  order  statistic. 

A  sample  mean. 

A  channel  in  the  multichannel  radar  return  xin). 

The  state  vector  process  in  the  state  space  model. 

In  Ozturk’s  algorithm,  tte  standardized  value  of  an  order  statistic. 

An  order  statistic  generated  in  Ozturk's  algorithm  by  Monte  Carlo  simulation 
of  the  null  distribution. 

A  matrix  in  the  Canonical  Correlatioas  Algorithm. 

(1)  A  loglikelihood  statistic.  (2)  The  diagonal  matrix  of  eigenvalues  in  the 
Singular  Value  Decomposititm  of  £. 

The  ith  order  statistic  in  a  sample  of  loglikelihood  statistics. 

The  conrelation  matrix  at  lag  /;  notation  used  in  die  Canonical  Correlations 
Algorithm. 

An  estimate  of  the  correlatkm  matrix  at  lag  /. 

The  diagonal  matrix  of  eigenvalues  in  the  Singular  Value  Decmnpositicm  of 

The  (zeroth  lag)  correlaticm  matrix  ammig  the  state  space  vectors. 

(1)  The  covariance  matrix  for  ein ).  (2)  The  covariance  matrix  for  x  in  ). 

An  estimate  of  the  covarianw  matrix. 

The  covariance  matrix  for  the  driving  noise  process  in  an  AR  model  of  the 
clutter. 


£m  The  covariance  matrix  for  the  white  noise  process  win). 

Q  The  (zeroth  lag)  correlation  matrix  fcx  the  multichannel  innovations  process 

in  the  innovations  rqxesentation  of  the  state  space  model. 
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(1)  A  probat^ty  defining  the  tail  in  the  Extreme  Value  Theory.  (2)  The 
MgftificRnce  level  of  a  statistical  test.  (3)  The  locatitm  parameter  of  the 
dutribution  to  be  estimated  by  Ozturk’s  algorithm.  (4)  A  parameter  used  in 
generating  K  distributed  SIRPi. 

The  level  of  a  statistical  test 

The  state  vector  in  the  innovations  representation  of  the  state  ^>ace  model. 

An  esfimate  of  the  state  space  vector  in  the  innovatioos  rqxesentation  of  the 
state  ^>ace  model. 

(1)  The  scale  parameter  of  die  distribution  to  be  estimated  by  Ozturk’s 
algocidim.  (2)  A  parameter  of  the  Gamma  probability  distributi<xt. 

The  state  vector  after  a  basis  transformation  of  the  innovations  representation 
of  the  state  space  model. 

(1)  A  parameter  in  the  Gm^alized  Pareto  distribution.  (2)  A  confidence 
level. 

An  estimate  of  a  parameter  in  the  Generalized  Pareto  distribution. 

A  probability  weighted  moment  used  in  the  Extreme  Value  Theory. 

An  estimate  of  cq. 

A  probability  wdghted  mmnent  used  in  the  Extreme  Value  Theory. 

An  estimate  of  £(. 

A  driving  noise  term,  uncorrelated  across  time  but  possibly  correlated  across 
channels. 

The  mean  over  time  samples  of  process  tin ). 

The  driving  noise  term  for  die  clutter. 

A  direshold  for  a  given  false  alarm  {vobability. 

An  estimate  of 

The  angle  with,  the  abscissa  of  linked  vectors  in  Ozturk’s  algorithm. 

A  loglikelihood  statistic. 

The  loglikelihood  statistic  defining  the  tail  in  die  Extreme  Value  Theory. 

The  tme-lag  temporal  correladoo  (matrix)  parameter  for  the  clutter. 

Eiqiected  values  of  certain  order  statistics  in  Oztuik’s  algorithm. 

Innovations,  uncorrelated  both  in  time  and  across  channels. 

Tlw  innovations,  uncorrelated  both  in  time  and  across  channels,  fw  the  null 
hypothesis  model. 

The  innovations,  uncmrelated  both  in  time  and  across  channels,  for  the 
alternative  hypothesis  model. 

A  zero-mean  normally  distributed  random  vector,  uncorrelated  both  in  time 
and  across  channels;  used  in  generating  Gaussian  white  nmse  win). 

The  jth  canonical  correlation,  sorted  in  decreasing  order,  in  the  Canonical 
Correlatimis  Algorithm. 
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A  parametBr  in  tl»  Gentfalized  Pareto  distribution. 

An  estimate  of  o. 

A  variance. 

An  estimate  of  a  channel  variance. 

A  parameter  in  a  rqMrameterizatioa  of  the  Generalized  Pareto  distributitm. 

An  estimate  of  t. 

The  Doppler  shift  for  the  clutter. 

Ttw  cumulative  probability  distribution  function  for  the  loglikelihood  statistic 
under  die  null  distribution. 

An  i^^oximation  of  the  cumulative  distribudcm  function  for  the  loglikelihood 
statistic  under  the  null  distributimi. 

The  cumulative  distribution  function  for  the  Generalized  Pareto  distribution. 

A  function  used  in  calculating  Ordered  Sandies  Least  Squares  estimates  of 
the  Generalized  Pareto  distribution. 

The  probability  density  functim  for  the  random  variable  S. 

The  probability  density  function  for  the  random  variable  V. 

The  probability  density  fimctitm  for  the  randran  variable  W. 

The  Gamma  function. 

The  Cumulative  Distributicm  Function  for  the  reference  distribution  in 
Ozturk’s  algorithm. 
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Appendix  B. 

A  Single-Channel  MiniSystem 

During  a  previous  effort  (Vienneau  93).  a  version  of  the  menu-based  system  was 
delivered  to  Dr.  Jwge  Rmneu.  Hiis  minisystem  contained  capabilites  for  synthesizing  single¬ 
channel  fvoceses  and  saving  them  to  a  file.  A  new  version  of  this  minisystem  was  produced 
under  th^  effort.  This  new  versitm  of  die  minisystem  is  distinguished  from  the  original  by  the 
addition  of  a  capability  to  calculate  a  quadratic  form  useful  in  assessing  the  performance  of 
SIRP  synthesis  routines. 

The  system  is  delivered  on  a  tar  ti^  containing  source  code  and  object  files  for  certain 
lil»‘ary  routines.  It  is  intended  to  be  read  on  a  Sun  computer.  Object  code  for  the  libraries  is 
jM'ovided  for  bodi  the  Sun  3  and  Sun  4  architectures.  Reading  the  tape  creates  a  subdirectory 
called  "schan"  under  the  user’s  current  directory.  This  subdirectory  contains  the  files 
described  in  Table  A-1. 

To  run  this  system,  there  must  be  an  environment  variable  called  ARCH  already  defined 
in  the  user’s  environment.  Currently  the  Sun  3  and  Sun  4  architectures  are  supported  and 
ARCH  must  be  set  to  one  of  these  two  values.  Also,  die  software  requires  access  to  the  Unix 
f77  and  make  programs. 

The  system  is  executed  by  first  changing  the  working  directmy  to  "schan.”  Inen  "schan 
[return]"  is  entered  at  the  Unix  prompt  to  invoke  the  Single  Charmel  Process  Synthesis  Menu. 
The  menu  contains  the  following  cations: 

•  Generate  K-distributed  SIRP 

•  Genmate  Gaussian  noise 

•  Generate  single  charmel  AR  process  with  SIRP  ot  Gaussian  driving  noise 

•  Calculate  the  SIRP  quadratic  form 

•  Calculate  histogram  data  for  the  SIRP  quadratic  form 

•  Convert  a  file  to  ASCII  in  a  format  suitable  for  Dr.  Romeu’s  software 

•  Convert  data  to  ASCII  file 

The  user  should  select  the  apprq)riate  menu  prompt  The  selected  programs  are  automatically 
cmnpiled  if  necessary. 
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Table  B*l:  Files  in  the  Minisystem 


File 

Description 

README 

Describes  how  to  install  and  run  the  system 

An  empty  subdirectory 

source 
souroe/cbta.f 
source/cgwiLf 
source/data2ascii.f 
source/qhistf 
source/qsirp.f 
source^scarsiip.f 
souice/siipmi 
source/Makefile 
source/sun4.make 
souroe/suii3  .make 

A  directory  of  Fortran  source  code. 

Converts  binary  data  files  to  ASCII 

Generates  single  dumnel  Gaussian  white  noise 
Ctmveits  binary  files  to  ASCII  in  Romeu  format 
Histogram  of  quadratic  form  for  SIRPs 

Quadratic  form  for  SIRPs 

Generates  single  channel  AR  processes 

Generates  single  channel  SIRP  noise 

Compiles  and  links  above  programs 

Called  firom  Makefile 

Called  from  Makefile 

schan 

A  c-shell  for  invoking  the  system 

hip 

hlp/cbta-hlp 

hlp/pgwn.hlp 

h4)/data2ascii.hlp 

h4)/scarsirpJi4) 

h^sirpmJilp 

A  directory  of  help  files  associated  with  tlw  main  programs 

suii3 

sun3/lib£Bl.a 

sunS/libmata 

sun3/libmc.a 

sunS/libmv.a 

sun3/libvec.a 

Sun  3  object  code  for  various  library  routines 

suii4 

sun4/Iibfsl.a 
sun4/libinata  | 
suii4/libmc.a 
sun4/Iibmv.a 
suii4/libvec.a 

Sun  4  object  code  for  various  lilx’ary  routines 

bin 

bin/sun3 

bin/sun4 

bin/exe 

bin/menul 

A  directory  to  cmttain  linked  executables  and  object  code 

A  directory  to  contain  cmnpiled  Sun  3  object  code 

A  directory  to  ctmtain  compiled  Sun  4  object  code 

A  directory  to  contain  linked  executi^les 

A  c-shell  for  die  main  menu 

Appendix  C. 

Generation  of  K-Distributed  SIRPs 

The  MSPSS  includes  an  option  for  generating  K-distributed  Spherically  Invariant 
Random  Processes  (SIRPs).  This  capability  was  added  during  a  previous  effort  (Vieimeau  93). 
The  al^xithm  specified  in  (Rangaswamy  91)  is  used  to  generate  a  SIRP.  This  algorithm  asks 
the  user  to  enter  the  two  parameters  a  and  b.  As  a  matter  of  fact,  SIRPs  can  be  generated 
with  knowledge  only  of  a.  The  parameter  b  is  totally  unneccesary.  This  appendix  provides  the 
madiematical  foundation  for  this  claim.  The  af^oaches  described  here  are  not  implemented  in 
die  l^PSS.  This  (qppendix  is  only  provided  to  document  these  findings. 

K-distributed  SIRPs  rely  on  the  generation  of  a  random  variable  5  with  the  following 
Probability  Density  Function  (PDF); 

5  >0.  (C-1) 

Rangaswamy  suggests  S  can  be  generated  based  on  the  following  theorem: 

Theorem:  Let  V'  be  a  random  variable  with  the  following  PDF; 

/,(v)-— v>0.  (C-2) 

Let  S  •Via,  where  a* «  Then  S  is  distributed  as  above  with  PDF  f,(s ). 

JW 

The  random  variable  V  can  be  generated  from  V  •  where  W  is  from  a  Chi  Squared 
distribution  with  2  a  d^rees  of  freedom. 

Alternately,  S  can  be  generated  based  on  either  of  the  following  two  theorems; 

Theorem:  Let  be  from  a  Chi  Square  distribution  with  2  a  degrees  of  freedom.  Let 
Then  5  is  distributed  as  above  with  PDF  f,(5 ). 


Theorem:  Let  be  from  a  Gamma  distribution  with  parameters  a  and  P.  where  p  >  l/o. 
That  is,  W  has  the  following  PDF; 

x>0.  (C-3) 

Let  5  ■  VvF  Then  S  is  distributed  as  above  with  PDF  /,(5 ). 
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Appendix  D. 

One-Pass  Algorithms  for  Calculating  Moments  of  Distributions 

Many  programs  in  the  Multichannel  Signal  Processing  Simulation  System  provide  the 
user  with  the  capability  to  estimate  certain  parameters  for  each  of  many  realizations  of  a 
stochastic  process.  Estimated  parameters  include  Autwegressive  coefficients,  means,  the 
covariance  matrix,  and  correlation  functions.  Often  the  program  reports  the  mean  and  variance 
of  these  estimates  calc.^dated  over  all  realizations  of  the  stochastic  process. 


D.l  Means  and  Variances 


The  structure  of  the  MSPSS  makes  it  convenient  to  calculate  these  means  and  variances 
from  one  pass  over  the  stochastic  process  realizations.  Since  the  variance  is  defined  in  terms 
of  the  mean,  calculating  the  variance  from  one  pass  through  the  data  requires  some  care.  The 
so-called  "calculator"  formula,  found  in  many  statistics  textbooks,  should  not  be  used: 


1 


^  (x,  -x)^ 


1 


in  -  1) 


EV- 

1  - 1 


(n  -1) 


xV 


(D-1) 


Equation  D-1  yields  a  one-pass  algorithm  for  calculating  the  variance,  but  this  algorithm  is 
numerically  unstable.  The  two  terms  in  the  difference  can  be  of  the  same  order  of  magnitude, 
leaving  the  variance  to  be  detennined  by  round-off  error. 

The  approach  for  calculating  the  variance  was  obtained  from  (Chan  83)  and  is  based  on 
the  two  following  theorems: 


Theorem  D-1:  Suppose 

A/t -A/».,  +  |(x* -A/*.,)  (D-2) 

and 


A/o-0.  (D-3) 

Then 

(D-4) 

In  other  words,  as  defined  by  the  above  difference  equation,  is  the  sample  mean  of  x,,  xj, 

...,  x*. 


Theorem  D-2:  Let  be  as  above.  Suppose 


(k  -  1 ) 5*2  -  (k  -  2)5*2.  ,  +  (X*  -  Af*  )2. 

(D-5) 

fw  it  =  2,  3,  4,  ...  Suppose 

5,2  -0 

(D-6) 

Then 
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(P■^) 

for  k  s  2, 3. 4. ...  In  other  words.  St  is  the  sample  variance  of  X|.  X2 . ji*  ■ 

A  c<»nputer  program  based  on  these  thewems  would  update  Mt  and  St  from  a  single  pass  of 
the  date.  Mt  must  be  updated  before  St  since  Mt,  not  Mt.i,  ^)pears  in  Equation  D-S.  On 
completiai.  the  final  values  of  Mt  and  St  will  be  the  desired  means  and  variances. 

For  completeness,  proofs  of  these  two  theorems  are  provided  below.  A  discussion  of 
their  numeric^  properties  is  provided  in  the  referenced  paper. 

Proof  of  Theorem  D-1  (by  induction): 


Hist,  show  the  theorem  holds  for  k  •  1.  By  hypothesis. 
Wi-Afo+ j(x,-Afo) 

Afi  ■  0  +  l(xi  -  0)  -Xj 
I  1 

Ml  -T  I 

*  i  -  J 


(D-8) 

(D-9) 

(D-10) 


Second,  assume  the  thenem  hteds  for  it  -  1: 

Ji  .  (D-11) 

*  ”  ^  i  ml 


Hnally.  prove  the  themem  holds  fx  k : 
Mt  “  A/*  _  1  +  —  (xj  -  A/*  _  I ) 

Af*- j(*Af*.,+x*-Af*.,) 


Aft  -j[(A-l)Af*.,+Xt) 


Aft  - 


k  ~  I 


Af* 


k 


'k  - 1 

z  •*- 


+  Xt 


Aft 


*  I  - 1 


But  this  is  what  was  to  be  showa 


(D-12) 

(D-13) 

(D-14) 

(D-15) 


(D-16; 


(D-17) 


78 


Proof  of  TlMorem  D-2  (by  ioducdon): 


Titsx,  show  the  tbeoiem  holds  for  1:  ■  2.  By  hypothesis, 
(2  -  1)S|  -  (2  -  2)S|.  I  +  (jcj  -  Af j)* 

1S|  -OSi  +  j[jr2-  + 


Si  -2 


f  'll  C  12 

X,-Jt2  . 

S2  -  -3—  +  — j— 


2  (  12 

JCi  -I- 


X,+X2  JC,+X2 

- —\  - — 

Si  m  (Xi  -  A/j)*  +  (Xj  - 

£(■*.  - Af 2)' 

^  <  ■  1 


(D-18) 


(D-19) 


(D-20) 


(D-21) 


(D-22) 


(D-23) 


(D-24) 


Secxmd.  assame  the  theorem  holds  for  A  -  l: 

Heooe. 


(D-25) 


(*-2)S»^,  -  £(x,  -Af*..)* 

I  - 1 


(D-26) 


Rnally,  prove  the  theorem  holds  for  k  : 

£  (X,  -Af»)2-*£  (Xj  -Aft)2  +  (x»  -Af»)* 
i  >  1  i  >  1 

£  ( X,  -  Aft )»  -  ( X,  -  Af*  . ,  -  I  ( X*  -  Af*  . , )  ]* 

i  ml  i  ml  < 


(D-27) 


(D-28) 
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ZiXi  -A/*)*-  Z  u Z  If  -Wi.i) 

<•1  i-l 

* -1 f  1  v 


(*  -1) 


L  u 

i  -  1 


L  (JC. 

j  -  1 


-  A/*  )*  -  (*  -  2)5*^  1  -  jU  -  A#*  . ,)  SU 

^  I  ■  I 

+  -Aft-i)Aft-i-»-  '*  ~/\xt  -A/t.i)" 

-  Af*  )*  -  (t  -  2)S^ ,  -  (X*  -  Af,  .  , )  £’x, 


+  (It  -1)  ^■'*  (^  -V)V  '►A#t.,)2 

-Af*)*-(*-2)5t^,  -  -Aft.,)A/*., 

+  U  '  A/t . , )  A/*  . ,  +  (X*  -  4/t  . ,  )* 

-  (It  17)^2  -•**)- (A Af*.,- Aft _,)f 

.A#*)*-a-2)5*»., 

•Aft)*-(t-2)5t*.,  +-^j^U-Aft)* 


(A-1) 


Lu 

1-1 


L  (•>?? 

- 1 


Z  (^i 

m  1 


U  -  Aft  - ,  )2  -  ( xt  -  Aft . ,  )* 


(D-35) 


£(jr. -Aft  )*-(!: -2  )Sr.,  +  77-^7  (x*  -  Af*  )* 
i  >  1  (,«  -  1  j 

Tberefore. 


( A  -  1  )St*  -  ( A  -  2)St*. ,  +  )*  (D-36) 

a-l)5**-  i(x, -Af*)*  (D-37) 

i  -  1 


Or 


tf-flTTyfu 

which  was  to  be  shown. 


(D-38) 


DJt  Correlation  Coefficients 

A  one-pass  algtnithm  for  calculating  the  correlation  coefficient  was  needed  in 
implementing  Ozturk’s  algorithm.  The  algorithm  implemented  relies  on  the  following 
theorem: 


Theorem  D-3:  Suppose 

53tyi-AfJ -AfJ-O,  (D-39) 

and 

Aff  -  Af/. ,  +  i-  (ai  -  Aft'., ) .  (EMO) 

Aft  ■Af/.i  •|•(3rt -Aft_, ),  (D-41) 

for  it  s  1,  2.  3. ...  Suppose 

sxy*  -  an  - ,  +  (D-42) 

for  A  s  2.  3. 4, ...  Then, 

ay*  -  £  (Xi  -  Aft'Xy,  -  Af/)  (I>43) 

i  ■  1 


Given  two  random  samples,  x,.  xj, ....  x,  and  3>,,  y2.  J'n.  the  customary  sample  estimate 
of  the  correlatitm  coefficient  is  given  by  Equatitm  dUO; 


5Xy. 


(«  -  DS'S: 

A  proof  of  Theorem  D-3  follows. 


(D-44) 
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Proof  of  Theorem  D<3  (by  induction): 


From  Theorem  D-1.  we  know 

*  i  •  i 

t 

*  i  -  1 

Fust,  show  the  theorem  bolds  for  it  •  2.  By  hypothesis. 

5X72  -  SXy,  +  (Jtj  -  Aff  Xyj  -  ) 

5X72  - 5X7,  +  Y  [xj  -  i  (j:,  +  Jt2)](y2  -  y  (y I  +  y2)l 

-[ -j  (*1  + -  •Jfalt  Y  (yi  +  yj)  -  yil  +  Ijt2  -  j  (jfi  +  *2)1  [>2  “  y  (y  1  +  72)1 


5X72- 


2x,  -  x,  +  X2  -  2x2 


2yi-yi  +  y2-2y2 


+  (x2-A/f  )(y2-A/5) 


5X7, 


»  » 

X,+X2 

71+72 

■*1  “  — — 

2 

2 

+  (jf2-Af?)(y2-Af2) 


SX72 - (X,  -  A/f  )(y,  - ilf|)  +  (X2 - Aff  Xyj - AfJ ) 
3*^2- EU -AffXy, -AfJ) 


i  -1 


(D-45) 

(D-46) 

(D-47) 

(D-48) 

(D-49) 

(D-50) 

(D-Sl) 

(D-52) 

(D-53) 


Seamd,  assume  the  theorem  bdds  for  A  -  1: 
SX7t . ,  -  e‘  (Xi  -  Af?, ,  Xy,  -  Af/. ,  ) 

I  - 1 


(D-54) 


Finally,  prove  the  theorem  bolds  for  A; 

£  (X,  -  Af/Xy,  -  A#/)  -  ‘e  -  A#/)  +  (X*  -  Mhiyt  -  Ml)  (0-55) 

i  ■  1  i  •  1 


E  -  Ml){yi  -  Ml)  -  *E  -  Ml. ,  -  I  (Xt  -  Ml. ,  )][y,  -  Ml. ,  -  j  (y*  -  Ml. ,  )] 

1-1  1.1  *  * 


+  (x*  -Ml)(yt-Ml) 


(D-56) 


SZ 


E  iXi  -  MhiVi  -  -  *E  ( -  wf- ,  Hy'i  -  ML , ) 

1-1  I « 1 

-7U*  -A/**-i)  e‘u  -  a//.,) 

*  1-1 

-yCy* -i^/.i)*!  U -w**.,) 

^  i  ■  1 

+  (Jt*  -  W/.  1  Hyt  -  ML  I ) 

+  (jtt  -ML(yt-ML 

The  inductive  hypothesis  is  used  at  this  step. 

£  (jti  -  MLiyi  -  A//)  -  sxy* -  7  (jf*  -  ML  1 )[  £*y,  -  (*  -  1  )Af^^ ,  i 

i  -  I  *•  1-1 

-7(3'*-  ML  I  )[*£  X.  -  (t  -  1  )Mf.  ,  ] 

^  i  -  1 

+ » ><>*  - •  > 


£  (x,  -  A/**)(y,  -  Af/)  -  5X7*  _ ,  -  I  (x*  -  Af**.  i ) ((*  -  1  )Af*''. ,  -  (it  -  1 ) Afi'. ,  ] 

i  -  I  * 

-  {(7* -Afi'.,  )[(*-! )A#**.,  -(jfc-DW*.,] 

+  (X*  -  Aff. ,  )(y*  -  A/r , ) 

+  (x*  -MLiyk  -  ML 

£  (X,  -  MLiyi  -  ML  -  sxy*  - 1  +  (x*  -  Aff . ,  )(y*  -  ml  1 ) 

1-1  * 

+  (xt  -MLiyt-ML 


(D-57) 


(D-58) 


(D-59) 


(D-60) 
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(D-61) 


f  (Xi  -  -  O  -  SXTt .  i  +  (x*  -  ,)()-*-  A//. , ) 


-Mh(yt  -Ml) 

--fJ^(Xk-Ml)(yt-Ml) 

L  U  -  Ml){yi  -  M/)  -  SXy»  . ,  +  ■■  *-  .(xt  -  Ml)iy,  -  Ml) 

im\  i*  •  1/ 


+  (**  -  Ml. ,  )(>*  -  Ml. , ) 

-  '  W/.i  -j(Xi -Ml.,] 


[yt-Ml.,-j(yt~Ml.,] 


-Ml)(y,  -Ml)mSXy,.,-h 


■  - 1 


j^-^(xt-Ml)(yt  -Ml) 


+  . ,  )(y*  -  Ml. ,  ) 


1 


(*-l) 


k  ^  k  *•* 


t  yt'  t  "*-* 


i  U  -  Ml)(yi  -  Ml)  -  sxy*  . ,  +  Ji-^ix,  -  Ml)(y,  -  Ml) 


i  m  1 


+  (xt  -  Ml.  ,){y,-  Ml. ,  ) 


JJ-^^^^i*k-Ml.,]^^^[y,-Ml.,] 


Unis. 

z  (Xi  -  MDiy,  -  Ml)  -  sxy*  . ,  +  (X*  -  Ml)(y,  -  Ml), 

i  ml  X  -  1 

which  was  to  be  shown. 


(D-62) 


(I>^3) 


(D-64) 

(D-65) 
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Appendix  £. 

RANDOM  VARIATE  GENERATION 

Ozturk’s  algorithm  relies  on  the  simulation  of  random  variates  from  certain  probability 
distributioos.  Currently,  15  distributions  are  implemented  in  the  software  inq}lementation  of 
Ozturk’s  algorithm,  lliis  ^)pendix  defines  these  distributions  and  explains  how  variates  are 
generated  from  them. 


E.1  Uniform  Variates 

Variates  from  aU  other  distributi(»s  are  generated  based  on  transformations  of  variates 
from  a  uniform  distribution.  For  example,  a  common  technique  for  generating  a  variate  frmn  a 
given  distribution  is  to  invert  the  Cumulative  Distribution  Function  (CDF).  Then  a  random 
variate  X  can  be  generated  from  the  desired  distribution  by  Equation  E-1: 

(E-1) 

where  U  is  uniform  on  (0,  1)  and  F  is  the  CH^F  of  the  desired  distribution. 

The  generation  of  random  variates  fr*om  a  uniform  distribution  is  a  problem  commonly 
treated  in  the  literature.  The  random  number  generator  used  in  the  MSI^S  is  based  on  an 
"universal"  random  number  generator  algorithm  develq)ed  by  Geor^  Marsaglia  of  the 
Supercomputer  Computations  Research  Institute  (SCRI)  at  Florida  State  University.  This 
algorithm  has  passed  stringent  tests  fw  randomness  and  indq>endence.  The  generator  has  an 
extremely  long  period,  about  2'^,  and,  fcsr  the  default  seed  values,  it  generates  exactly  the 
same  sequence  of  24-bit  random  numbers  in  a  variety  of  languaga  and  machii^.  (Harmon  88) 


EL2  Beta  Variates 

A  method  frx  generating  variates  frmn  Beta  distributions  was  already  embodied  at  the 
beginning  of  this  in  the  Gamma  generate  for  the  MSPSS.  This  mefiiod,  which  differs  frmn 
that  used  by  Ozturk,  is  described  here. 

The  Probability  Density  Fimction  (PDF)  of  a  Beta  distribution  is  given  by  Equation  E-2: 


The  mean  of  a  Beta  distribution  is 


a  + 1 


',  and  the  variance  is 


(a  -t-lXh  +1) 


a+b+2  (a +1> +2)*(a +6 +3)‘ 

The  method  used  for  transforming  uniform  variates  into  Beta  variates  is  based  on  the 
following  theoran; 

Theorem:  Let  U i  and  C/2  be  independent  identically  distributed  random  variables  from 
a  uniform  distribution  on  (0,  1].  Let 


r,  -  u}'“ . 

(E-3) 

(E-4) 

If  Ki  -t-  y2  is  less  than  or  equal  to  unity,  set  Y  to 

(E-5) 

Then  Y  is  from  a  standard  Beta  distribution  with  parameters  a  and  b. 
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1 


n 

'  N  • 

2  J 

EJ  Cauchy  Variates 

The  PDF  of  the  Cauchy  distributicm  is 

The  Cumulative  Distributitm  Fuoctitm  (CDF)  is 

2  ft 

A  variate  can  be  generated  frtmi  this  distribution  as 

X  m  Tan 

where  U  is  uniform  on  (0,  1). 

E.4  Exponential  Variates 

The  PDF  for  the  Exptmential  distribution  is 

X  iO. 

The  CDF  is 

Fix)ml^e-'^.  X  iO. 

A  variate  from  this  distributitm  can  be  generated  as 

X  —iiLn(U), 

where  U  is  uniform  tm  (0,  1]. 

E^  Extreme  Value  Variates 


(E-6) 


(E-7) 


(E-8) 


(E-9) 

(E-10) 

(E-11) 


The  PDF  for  the  Type  I  Extreme  Value  distribution  for  minimum  elements,  also  called 
the  Gumbel  distribution,  is 


1 


where  o  is  the  scale  parameter.  The  CDF  is 


F(x)-  l-e-'. 


(E-12) 


(E-13) 


(Dr.  Ozturk  code  is  mistaken  in  the  calculation  of  the  Extreme  Value  CDF  in  the  routine 
THETA  for  calculating  the  angles  between  the  linked  vectors  and  the  abscissa.)  A  random 
variate  can  be  generated  fr<«i  this  distribution  as 


X  maLn[-LniU)]. 


(E-14) 


where  U  is  uniform  on  (0,  1). 


E.6  Gamma  Variates 

The  capacity  to  genorate  Gamma  variates  already  existed  in  the  MSPSS  at  the  beginning 
of  this  effort.  The  algoridim  used  is  taken  from  (Fishman  73).  It  tends  to  create  overflow  and 
underflow  {problems  for  a  <  0.25  or  a  >  100.0.  Warning  messages  alert  the  user  to  these 


problems.  The  algorithm  tends  to  beoHne  much  slower  as  the  number  of  degrees  of  freedom 
increases.  This  algorithm  differs  from  the  one  embedded  in  Ozturk’s  code. 

The  PDF  for  a  Gamma  distributitm  is 

where  a  is  the  shape  paranoettf  and  ^  is  the  scale  parameter.  The  mean  of  the  Gamma 
distributioa  is  and  tte  variance  is  af)^. 

The  algoridim  first  generates  a  random  variable  Y  with  a  Gamma  distribution  with  a 
scale  parameter  equal  to  unity.  Once  Y  is  generated,  a  variable  X  is  set  to  fiY.  X  that  has  the 
desired  sluqje  and  scale  parameters.  This  procedure  relies  on  the  following  theorem: 

Theorem:  Let  X  be  a  random  variable  from  a  Gamma  distribution  widi  an  arbitrary 
shape  parameter  a  but  a  scale  parameter  equal  to  one.  Let 

X-pK.  (E-16) 

Then  X  has  a  Gamma  distribution  with  shiqje  parameter  a  and  scale  parameter  p. 

The  random  variate  Y  is  generated  by  particming  a  into  an  integer  part  and  a  real  part. 
Gamma  variates  are  generated  for  each  of  these  parts,  and  their  sum  is  the  desired  Gamma 
variate: 

Theorem  Let  V  be  from  a  Gamma  distribution  with  a  shape  parameter  k  •  LaJ  and  an 
unit  scale  parameter.  Let  W  be  fr<»n  a  Gamma  distribution  with  a  shape  parameter 
a  ~-k  and  an  unit  scale  parameter.  Let 

YmV^W  (E-17) 

Then  Y  is  from  a  Gamma  distribution  with  shape  parameter  a  and  unit  scale  parameter. 

A  Gamma-distributed  random  variate  with  an  integer  shape  parameter  and  unit  scale 
parameter  is  generated  as  the  sum  of  exponential  random  variates: 

Theorem  Let  At,  i  =  1,  2,  ...,  k  be  ind^)endent  and  identically  distributed  randtmi 
variables  from  an  e^qxmential  distribution  witii  unit  means.  Let 

■  Aj  +  Aj  +  •  ■  •  +•  Aj  .  (E-18) 

Then  V  is  Gamma  distributed  with  a  sh^  parameter  of  k  and  an  unit  scale  parameter. 

A  Gamma-distributed  random  variate  with  a  real  shape  parameter  between  zero  and  one 
is  generated  as  the  product  of  a  Beta-distributed  random  variate  atxl  an  erqxmentially- 
distributed  random  variate: 

Theorem  Let  Z  be  from  a  Beta  distribution  with  parameters  a  =  a  and  b  s  1  -  a.  Let 
A  be  e3q)onentially  distributed  with  unit  mean.  Let 

W'-ZA.  (E-19) 

Then  W  is  Gamma  distributed  with  shape  parameter  a  and  unit  scale  parameter. 

E.7  Gumbel  (Type  II)  Variates 

The  PDF  fm*  a  Gumbel  (Type  11)  Extreme  Value  distribution  is 
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fix) 


-1 


-f-i  _ 

M 

X 

a 

e 

a 

-» 


.  X  >0. 


where  a  is  the  scale  parameter  and  p  is  the  shape  parametm*.  The  CDF  is 


F(x)»e 


.  X  >0. 


A  variate  can  be  gmerated  from  this  distribution  as 

-1 


Xma 

where  U  is  uniform  on  (0,  1). 

EJ  Johnson  Variates 

The  PDF  of  a  Johnson  Su  distribution  is 


LniU) 


(E-20) 


(E-21) 


(E-22) 


f(x)~ 


^2nxH(y  - 


Ejq> 


1 

2^1 


Y-6  +  TiLrt{  ^-Lv  -  e  +  V(y  -  e)^  +  ) 


(E-23) 


e  is  the  location  parameter,  A.  is  the  scale  parameter,  and  y  and  r\  are  shape  parameters.  A 
random  variable  from  a  Johnson  Su  distributimt  can  be  generated  by  transforming  a  standard 
normally  distributed  variate: 


X  «E4-Xsiiih| 


Lzx 


(E-24) 


where  Z  is  standard  Gaussian. 

E.9  K  Distributed  Variates 

The  ability  to  generate  random  variables  from  a  K  distribution  existed  in  the  MSPSS 
when  this  project  began,  .^ipendix  C  describes  how  this  capability  is  inqilemented.  The  PDF 
of  a  K  dishibution  iroludes  a  shiyie  parameto:  a. 

E.10  Laplace  Variates 

The  PDF  for  die  Laplace  distribution  is 


The  mean  of  this  distribution  is  zero,  and  the  variance  is  two.  The  CDF  is 

X  <0 


F(x)- 

A  variate  from  die  Liqilace  distribution  can  be  generated  by 


1  xiO 

2  • 


(E-25) 


(E-26) 


Ln(2U), 

-Ln[2(l-U)]. 


(E-27) 


whoe  U  is  uniform  on  (0,  1). 
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E.11  Logistic  Variates 

The  PDF  for  Ae  Logistic  distnbuti<»>  is 


The  CDF  is 


F(x)m 


1 


i  4- e~ 


A  vsriate  can  be  generated  from  this  distribution  as 

1-1/ 


X  —Ln 


U 


where  U  is  unifotn  on  (0,  1). 


(E-28) 

(E-29) 

(E-30) 


E.12  Lognormal  Variates 

A  random  variable  is  from  a  Lognormal  distribution  if  its  natural  logarithm  is  normally 
distributed.  The  PDF  of  a  lognormal  distribution  is 


fix) 


_1 

V2nax 


X  >0. 


fE-31) 


4  and  o  are  the  mean  and  standard  deviation,  respectively,  of  the  normal  transformation,  not 
the  mean  and  standard  deviation  of  the  lognormal  distribution. 

Lognormal  variates  are  generated  by  transforming  normal  variates: 

Xme^,  (E-32) 


where  Z  is  normally  distributed  with  mean  4  and  variance  o^.  (The  generation  of  lognormal 
variates  is  not  implmnented  coneedy  in  Dr.  Oztuik’s  code.) 


EL13  Normal  Variates 


The  capability  to  generate  normal  variates  already  existed  in  the  MSPSS  at  the  beginning 
of  this  dfort.  The  algorithm  described  hoe  is  based  on  (Press  86).  The  PDF  for  die  Normal 
distribution  is 


fix) 


1 

'l2na 


e  . 


wtore  4  is  the  mean  and  is  die  variance. 


(E-33) 


In  generating  normal  variates,  it  is  sufficient  to  design  an  algorithm  to  generate,  zero- 
mean,  unit-variance  normal  variates.  Normal  variates  with  arbitrary  means  and  variance  can 
dum  be  generated  by  means  of  the  following  transformation: 

X  mn  +  aZ.  (E-34) 


where  Z  is  normal  with  mean  zero  and  unit  variance.  The  resulting  variate  X  is  normally 
disiributed  with  mean  4  and  variance  o^. 

The  algoridim  used  to  ^nerate  standard  normal  variates  is  as  follows.  A  pair  of  random 
variates  (Ki,  Vj)  are  gmierated  where  dre  pair  is  uniformly  distributed  on  the  unit  circle.  Then, 
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die  random  variables  Z|  and  Zj  given  by 

Z.-v.-y/n^^  (E-35) 

(E.36) 

are  stochaatically  indqiendent  and  normally  distributed  with  mean  zero  and  unit  variance, 
where  R  is 

«  -  V?  ♦  Vi  .  (E-37) 


E.14  Pareto  Variates 

The  PDF  fw  the  Pareto  distribution  is 


The  mean  is  6162/(62-1). 
distributi<m  is 


/(j:)-6,'*^e2j'‘*''’‘\  Jti6,. 

The  variance  is  6f62/[(62- 1)*(62- 2)1. 


Fix) 


X  i6,. 


A  variate  frmn  the  Pareto  distribution  can  be  generated  as  follows: 

X  -6, (1/1/ )*'“*. 


(E-38) 

The  CDF  for  this 

(E-39) 

(E-40) 


where  U  is  uniform  (Hi  (0,  1], 


E.15  WeibuU  Variates 

The  PDF  for  the  WeibuU  distributi(»  is 

/(x)--^x«»-‘r-^"“^  xiO.  (E-41) 

whoe  a  is  the  scale  parameter  and  b  is  the  shape  parameter.  The  CDF  is 

f  (X )  - 1  -  X  2  0.  (E-42) 

The  mean  is  ar(l  1/6)  and  the  variance  is  o^[r(l  •►Z/P)-!^!!  1/6)1.  A  WeibuU  variate 
can  be  generated  as 

X -al-LB(^/)]*'^  (E-43) 

where  U  is  uniform  on  (0,  1}. 
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Appendix  F. 

CONFIDENCE  REGIONS 

Implementing  Ozturk’s  goodness-of-fit  test  requires  the  calculation  of  confidence  regions 
for  Ozturk’s  statistic.  This  appendix  outlines  the  theory  behind  these  calculations. 

F.1  Bivariate  Normal  Distribution 

Consider  the  problem  of  determining  a  confidence  ellipse  for  a  pair  of  raiulom  variables 
from  a  bivariate  normal  distribution.  Let  Xi  and  have  the  joint  Probability  Density  Function 

/  (xi.xj): 


(F-1) 


2RO1O2V1  - 


2 

^  1a 

^2-14 

11 

“1 

“2 

“2 

Jj 

This  is  the  PDF  of  a  bivariate  normal  distribution  where  X/  has  mean  and  variance  a,’,  i  = 
1,  2.  The  cmrelation  coefficient  for  Xi  and  Xj  is  p.  The  problem  is  to  determine  a  confidence 
region  C  such  that  the  probability  of  Xt  and  X2  being  in  C  is  100  y%. 

Consider  the  transfcxmation  defined  in  Equations  F-2  and  F-3: 

Z,-«,(X,.X2)-^i;^  (F-2) 


Z2  ■  «2(^1-^2)  ■  - Z -  “P  - Z - 

Vl  -  p* 

~  U  V  J  V  J  J 

The  inverse  transfomation  is  given  by  Equations  F<4  and  F-S: 

“  wi(Zi.Z2) "  Pi  +  OiZi 

X2  ■  w2(Zi.Z2)  “  P2  +  <r2pZi  +  02VI  -  p*Z2 
The  Jacobian  of  fins  transformation  is 


■  0x02^  1  -  • 


(F-2) 

(F-3) 


(F-4) 

(F-5) 


dwi 

dwi 

dzi 

dzj 

dW2 

dw2 

dzi 

dZ2 

(F-6) 


^  8  (21.Z2)  be  the  joint  PDF  fcx*  Z|  and  Z2.  From  a  theorem  in  analysis  relating  to  a  change 
of  variables  in  int^radon  (Hogg  78).  the  joint  PDF  of  Z,  and  Z2  is  related  to  the  joint  PDF 
for  X]  and  X2  by: 


g  (21.22)  -  I/W  1/  (»Vi(Zi.22).M'2(2,.22)). 


(F-7) 


*  (Zi.Z2)- 


1 


(F-8) 


So  the  transformadon  has  converted  the  original  random  variables  into  stochasdcally 
indq)endent  standard  nwmal  variates. 

A  confidence  region  in  the  transformed  space  is  a  circle: 
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a  (C  )  -  1  (zi.rj)  1  z?  +  <  c* ) , 

(F-9) 

wlmre  the  radius  c  is  chosen  so  as  to  ensure  the  {vobability  of  membership  in  this  circle  is  the 
desired  value  y.  This  probability  can  be  found  as  follows: 

Pr(zf  4-z}  <C^)m  j  ^  f  (Zi,Zz)dZi<b2 

(F-10) 

2n  e  _r^ 

Fr(z? +Zj  <c*)- 1  ^rdrdQ 

(F-11) 

J'r  (zf +2|  <  c*)  -  1  -  « 

(F-12) 

Set  this  ixobability  to  y,  and  solve  for  c : 

(F-13) 

c  -V-21n(l-y). 

The  confidence  region  for  Zi  and  7]  can  be  transformed  back  into  the  original  space.  Thus,  a 
100  ctmfidence  ellipse  for  and  Xj  is 

C  -  \  (JCi.JCj)  1  <  -21n(  1  -  y) ) .  (F-14) 

An  a{^>roximate  confidence  ell4)se  is  formed  by  using  sample  estimates  for  the  means, 
variances,  and  correlation  coefficients  in  Equation  F-14. 


F2  The  Johnson  System  of  Transformations 

Ozturk’s  statistic  need  not  be  from  a  bivariate  normal  distribution.  The  problem  then 
arises  of  approximating  a  confidence  region  for  other  distributions,  even  if  the  distribution  is 
unknown.  The  Jtfiinaon  system  of  transformations  addresses  this  problem  (Shah  93). 

A  distiibutitm  can  be  completely  characterized  by  its  moments,  when  they  exist  This 
suggests  that  any  distribution  can  be  dosely  apix’oximaled  by  a  distribution  matching  the  first 
few  mmnents.  The  Jdmson  system  of  transformatinis  ajyroximates  the  first  four  moments  of 
a  randmn  variable,  known  as  the  mean,  variance,  skewness,  and  kurtosis.  The  Johnson  family 
of  distiibutimis  is  defined  by  Equafitm  F-15: 


R  "y*r\fiiG;k,e),  i -1,2,3, 


(F-15) 


where  R  is  firom  a  standard  normal  distribution,  e  is  a  locatioo  parameter.  A.  is  a  scale 
parameter,  and  r  and  are  sluqpe  parameters. 

The  Jdmson  Sa  distribution  is  defined  by 


/,(G;X.e)-siiih-‘ 


G  -e 


The  Johnson  8$  distribution  is  defined  by 

G  -e 


/j(C;X.e)-lii 


^  +  e  —  G 


,  e<G  <e-»'A,. 


(F-16) 


(F-17) 


Finally,  the  Johnson  Si  distributitm  is  given  by 


/j(G;X.e)-ln 


G  »e 
X 


G  >e. 


(F-18) 


The  Johnscm  Si  distributitm  can  be  leparameterized  in  terms  of  y*.  where 
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(F-19) 


V* -Y-illn(e). 

The  transfomiaticm  for  the  Johnson  Sl  distribution  becomes 

i? -T  +tlIn(G -e).  G>e.  (F-20) 

The  Johnson  Si,  divides  the  skewness-kurtosis  space  into  two  parts.  The  Su  distribution 
occupies  one  of  the  resulting  regitMis,  and  the  distribution  occupies  the  other.  One  could 
which  member  of  the  Johnson  family  best  fits  a  random  sample  based  on  the  sample 
skewness  and  kurtosis.  Estimates  of  skewness  and  kurtosis,  howevo-,  are  highly  variable  for 
small  saiiq>les.  A  more  dqiendable  method  has  been  found  based  on  the  heaviness  of  the  tails 
as  compart  with  the  center. 

F^l  Fitting  the  Johnson  Distribution 

Let  G(i)  ^  G(2)  ^  ...  ^  G(,)  be  an  observed  ordered  sample.  The  problem  is  to  fit  G  with  a 
Johnstm  distribution.  Consider  any  given  positive  number  r  (0.775449  is  used  in  the  Ozturk 
algorithm),  r  is  used  to  divide  the  axis  for  the  normal  transformed  data  into  three  intervals, 
each  of  length  2r:  (-3r.-r),  (-r.r),  (r,3r).  Tltt  Johnson  family  of  transformatitxis  mtq>s  the 
observed  sample  into  percentile  ^timates  along  the  normal  axis:  Rn-,  i  Ri2)  ^  ^  R(n) 

^)I»'opriate  indices  for  approximating  the  endpoints  of  the  given  intervals  can  be  found  by 
solving  Equation  F-21: 

Pr(R<a)«^::^.  (F-21) 

n 

where 

a  m  -3r.  -f .  r .  3r .  (F-22) 

The  solution  is 

k,m[nPr(R  <a)  +  ll2}.  (F-23) 

where  [*]  denotes  the  nearest  integer.  Since  the  Johnstn  transformations  are  strictly 
monotonically  increasing,  the  corresponding  percentiles  in  the  untransformed  data  can  be 
estimated  as  gk_^,  gk_^,  gt^,  and  gt^^.  Parameters  that  characterize  the  tails  are  given  by 


Equations  F-24  through  F-26: 

fn  -&J, 

(F-24) 

(F-25) 

P  •gk,- 

(F-26) 

These  parameters  can  be  used  to  decide 
charactnizes  the  random  sample  as  follows: 

which  member  of  the  Johnson  family  best 

>  1  for  Johnson  Sy 

(F-27) 

■  1  for  Johnson  S- 
P 

(F-28) 
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<1  for  Johnsor 

Since  G  is  a  ontinuous  ^om  variable.  Equaticm  F-28  will  be 
Therefore,  in  practice,  one  $ccepta  the  data  is  from  a  Johnscm  St, 
side  is  in  a  small  region  around  unity. 

Parameter  estimates  for  the  Jtriinstm  Su  distributi(»  are: 


00811 

m  1 
—  +  — 

2 

P  P 

YaTjsiiih' 


,-i 


m 


y  p  p 


ipJsTT 

'  p  p 


y  p  p 


+  2 


Parameter  estimates  for  the  Johnson  St,  distributi<m  are: 


2r 


T 


e  • ' 


8r  +«. 


_£ 

2 


£L.l 


Parameter  estimates  fn*  the  Johnson  5b  distribution  are: 


r  _ 

cosh*' 

N 

1  +  ^ 
m 

. 

(F-29) 

true  with  iHX>bability  zero, 
distribution  if  the  left  hand 


(F-30) 


(F-31) 


(F-32) 


(F-33) 


(F-34) 


(F-35) 


(F-36) 


(F-37) 
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(F-38) 


yal^Sillh' 


,-l 


£_  £. 
/  m 

1  +  ^ 
ffl 

-4 

2 

ffl  / 
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1  +  f 

-2 

-4 

ffl 

2.L.I 
m  I 


gr  +«- 


2  2 


£-2. 
/  m 


£.£_1 
ffl  / 


(F-39) 


(F-40) 


Confidence  Regions 

A  confidence  region  for  Oztuilc’s  statistic  can  be  constructed  using  the  Johnson 
transformation.  The  M<xite  Carlo  data  for  (Un,  Vn)  is  used  to  determine  the  appropriate 
Johnstm  transformation.  Us  and  Vs  can  be  transformed  by  different  members  of  the  Johnstm 
family.  The  Monte  Carlo  data  is  then  transformed,  and  sample  means,  variances,  and  the 
correlatitm  coefficient  are  calculated  in  the  usual  manner  frmn  the  transformed  data.  A 
confidence  ell4>8e  is  constructed  in  the  transformed  Bivariate  Normal  space.  A  confidence 
region  in  the  original  space  is  then  formed  from  inverse  Johnstm  transformafions. 
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MISSION 

OF 

ROME  LABORA  TORY 


Mission.  The  mission  of  Rome  Laboratory  is  to  advance  the  sdence  and 
technologies  of  command,  control,  communications  and  intelligence  and  to 
transition  them  into  systems  to  meet  customer  needs.  To  achieve  this, 
Rome  Lab: 


a.  Conducts  vigorous  research,  development  and  test  programs  in  ail 
applicable  technologies; 

b.  Transitions  technology  to  currertt  and  future  systems  to  improve 
operational  capability,  readiness,  and  supportabiiity; 

c.  Provides  a  full  range  of  technical  support  to  Air  Fores  Materiel 
Command  product  centers  and  other  Air  Force  organizations; 

d.  Promotes  transfer  of  technology  to  the  private  sector; 

e.  Maintains  leading  edge  technological  expertise  in  the  areas  of 
surveillance,  communications,  command  and  control,  intelligence,  reliability 
science,  eiectro-magnetic  technology,  photonics,  signal  processing,  and 
computational  science. 


The  thrust  areas  of  technical  competence  include:  Surveillance, 
Communications,  Command  and  Control,  intelligence.  Signal  Processing, 
Computer  Science  and  Technology,  Bectromagnetic  Technology, 
Photonics  and  Reliability  Science. 


